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1.0  Introduction 

This  report  summarizes  the  results  of  an  assessment  of  the  accuracy  of 
response  of  six  numerical  step-by-step  procedures  used  in  computational  struc¬ 
tural  dynamics.  The  six  algorithms  used  in  this  study  are  representative  of  the 
different  types  of  numerical  procedures  used  to  compute  the  dynamic  structural 
response  to  a  time-dependent  loading  history.  The  time-dependent  loading  envi¬ 
sioned  in  this  research  is  that  of  the  motion  of  the  ground  below  a  discrete  struc¬ 
tural  model  and  is  expressed  in  terms  of  a  ground  acceleration  time-history.  The 
dynamic  structural  response  for  each  structural  model  used  in  this  study  is  charac¬ 
terized  by  the  computed  response  time-histories  of  accelerations,  velocities,  and 
displacements. 

All  structural  models  used  in  this  study  are  linear,  single-degree-of-freedom 
(SDOF)  systems.  The  natural  (undamped)  periods  Ta  of  these  SDOF  systems  are 
selected  based  on  consideration  of  the  important  modal  periods  of  hydraulic  struc¬ 
tures  such  as  gravity  dams,  arch  dams,  gravity  lock  walls,  U-ffame  locks,  and 
intake  towers.  The  forcing  functions  used  in  this  study  are  single  frequency  har¬ 
monics.  The  use  of  a  single  frequency  facilitated  the  evaluation  of  the  accuracy  of 
the  computed  responses  solved  for  at  regular  time  increments  during  ground 
motion. 

The  time  increments  A t  used  in  the  analyses  are  0.02,  0.01,  and  0.005  seconds. 
These  values  are  typical  of  the  At  used  in  discretizing  earthquake  acceleration 
time-histories  (e.g.,  Hudson  1979)  recorded  in  the  field  on  strong  motion  accelero- 
graphs  (shown  in  Figure  6.1.1  in  Chopra  1995). 
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1.1  Six  Methods  Used  to  Compute  Time  Domain 
Responses 

The  six  algorithms  included  in  this  study  are  the  Newmark  P  method  (with 
values  of  constants  y  and  P  corresponding  to  the  linear  acceleration  method),  the 
Wilson  0  method,  the  Central  Difference  Method,  the  4th  Order  Runge-Kutta 
method,  Duhamel's  integral  solved  in  a  piecewise  exact  fashion,  and  the  Piecewise 
Exact  Method  applied  directly.  All  of  these  algorithms  were  used  in  their 
discretized  forms  (i.e.,  the  loading  and  response  histories  were  divided  into  a 
sequence  of  time  intervals);  thus,  they  are  characterized  as  step-by-step 
procedures. 

The  six  algorithms  used  can  be  categorized  into  two  main  groups,  depending  on 
their  general  approach  to  satisfying  the  differential  equation  of  motion.  The  first 
group  includes  Duhamel's  integral  solved  in  a  piecewise  exact  fashion  and  the 
Piecewise  Exact  Method  applied  directly.  These  two  methods  formulate  exact 
solutions  to  the  equation  of  motion  for  assumed  forms  of  the  time-dependent  forc¬ 
ing  functions  (i.e.,  the  loading  is  approximated  by  a  series  of  straight  lines  between 
the  time-steps).  The  second  group  includes  Newmark  P  method,  Wilson  0 
method.  Central  Difference  Method,  and  4th  Order  Runge-Kutta  Method.  These 
methods  are  referred  to  as  numerical  methods  because  they  approximately  satisfy 
the  equation  of  motion  during  each  time-step  for  the  given  loading.  The  New¬ 
mark  P,  Wilson  0,  and  4th  Order  Runge-Kutta  methods  use  numerical  integration 
to  step  through  the  analysis  of  the  time  response  problem,  where  the  Central  Dif¬ 
ference  Method  uses  numerical  differentiation. 


1.2  Direct  Integration  Methods 

The  linear  acceleration  method,  Wilson's  0  method,  and  the  4th  Order  Runge- 
Kutta  method  are  examples  of  direct  integration  methods.  The  term  “direct” 
means  that  prior  to  numerical  integration,  there  is  no  transformation  of  the  equa¬ 
tions  into  a  different  form,  such  as  is  done  in  a  frequency  domain  analysis.  Inte¬ 
gration  methods  are  discrete  in  that  the  response  values  are  solved  for  at  regular 
increments  in  time  during  ground  motion,  which  are  separated  by  a  time  increment 
A/. 


Direct  integration  methods  are  based  on  two  concepts.  First,  the  equation  of 
motion  for  the  structural  model  is  satisfied  at  discrete  points  in  time  (i.e.,  t,t  +  At, 
t  +  2  At,  ...)  during  ground  motion.  Second,  the  forms  of  the  variation  in  displace¬ 
ment,  velocity,  and  acceleration  responses  within  each  time  interval  At  are 
assumed. 
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1.3  Difference  Between  implicit  and  Explicit 
Numerical  Methods 


Numerical  methods  such  as  direct  integration  methods  are  classified  as  either 
explicit  or  implicit  integration  methods.  Chopra  (1995),  Subbaraj  and  Dokainish 
(1989a,b),  Bathe  (1982),  and  Bathe  and  Wilson  (1976)  distinguish  between  the 
two  numerical  methods  as  follows.  The  explicit  integration  method  solves  for  the 
unknown  values  of  displacement  xt+ A/ ,  velocity  xt+ At ,  and  acceleration  x, + At ,  at 
each  new  time  t  +  At  using  the  equation  of  motion  for  the  structural  model  at 
time  t ,  with  the  unknown  values  for  xt,  xt,  and  xt  at  time  t  as  the  initial  conditions. 
The  implicit  integration  method  solves  for  the  unknown  values  of  xt+At ,  xt+ At , 
and  x*  j.  A  t  at  each  new  time  t  +  At  using  the  equation  of  motion  at  time  t  +  At.  For 
multiple-degree-of-freedom  (MDOF)  systems,  implicit  schemes  require  the  solu¬ 
tion  of  a  set  of  simultaneous  linear  equations,  whereas  explicit  schemes  involve 
the  solution  of  a  set  of  linear  equations,  each  of  which  involves  a  single  unknown. 
Thus,  the  explicit  integration  method  does  not  require  a  factorization  of  the 
coefficient  matrix  in  the  step-by-step  solution  of  the  equations  of  motion  for  the 
semidiscrete  MDOF  structural  system  model.  The  coefficient  matrix  is  a  com¬ 
bination  of  the  stiffness,  mass,  and  damping  matrices  of  the  MDOF  model. 

Along  with  others,  Subbaraj  and  Dokainish  (1989a,b)  noted  that  implicit  algo¬ 
rithms  are  most  effective  for  structural  dynamics  problems  (in  which  the  response 
is  controlled  by  a  relatively  small  number  of  low-frequency  modes),  while  explicit 
algorithms  are  very  efficient  for  wave  propagation  problems  (in  which  the  contri¬ 
bution  of  intermediate-  and  high-frequency  structural  modes  to  the  response  is 
important).  Accordingly,  of  the  two  types  of  numerical  methods,  implicit  algo¬ 
rithms  are  more  popular  in  earthquake  engineering  problems  because  of  the  larger 
time-step  that  may  be  used  in  the  analysis. 

However,  implicit  procedures  involve  considerable  computational  effort  at  each 
time-step  compared  with  explicit  methods  for  MDOF  semidiscrete  models  since 
the  coefficient  matrices  must  be  formulated,  stored,  and  manipulated  using  matrix 
solution  procedures.  Therefore,  in  blast  type  problems  where  a  small  time-step  is 
required  to  capture  the  structural  response  of  large-scale  models  involving  hun¬ 
dreds  to  thousands  of  degrees  of  freedom,  implicit  methods  are  computationally 
impractical,  and  explicit  methods  are  the  preferred  type  of  algorithm. 

Two  explicit  algorithms,  the  Central  Difference  Method  and  the  4th  Order 
Runge-Kutta  method,  are  included  in  this  study.  The  original  1959  linear  acceler¬ 
ation  method  version  of  the  Newmark  (3  family  of  numerical  methods  and  Wilson's 
0  method  are  classified  as  implicit  methods  (Newark  1959).  In  general,  implicit 
integration  methods  are  frequently  used  by  the  structural  dynamics/earthquake 
engineering  community  to  solve  for  the  response  of  semidiscrete  MDOF  structural 
models  to  earthquake  excitation. 
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1 .4  Questions  of  the  Accuracy  of  All  Six  Step-by- 
Step  Methods  and  the  Stability  of  Numerical 
Integration  and  Numerical  Differentiation 
Methods 


The  selection  of  the  size  of  the  time-step  At  to  be  used  in  the  step-by-step  cal¬ 
culation  of  the  dynamic  response  of  the  SDOF  (and  of  MDOF  semidiscrete  struc¬ 
tural  models)  is  restricted  by  stability  and/or  accuracy  considerations  for  the 
six  algorithms  included  in  this  study.  The  primary  requirement  of  a  numerical 
algorithm  is  that  the  computed  response  converge  to  the  exact  response  as  At  -  0 
(Hughes  1987).  However,  the  number  of  computations  increases  as  the  time-step 
At  is  made  smaller  in  a  dynamic  analysis,  an  important  issue  for  response  analysis 
of  semidiscrete  MDOF  structural  system  models. 

In  addition  to  accuracy  considerations,  stability  requirements  must  also  be  con¬ 
sidered  during  the  selection  of  the  time-step  At  to  be  used  in  a  step-by-step 
response  analysis  either  by  the  three  numerical  integration  methods  or  by  the 
numerical  differentiation  method.  Stability  criterion  is  expressed  in  terms  of  a 
maximum  allowable  size  for  the  time-step,  A tcritical.  The  value  for  A tcriticaX  differs 
among  the  four  numerical  algorithms. 

No  stability  criteria  (expressed  in  terms  of  a  limiting  time-step  value)  are 
needed  for  Duhamel's  integral  solved  in  a  piecewise  exact  fashion  and  the  Piece- 
wise  Exact  Method  applied  directly.  This  is  because  these  two  methods  formulate 
exact  solutions  to  the  equation  of  motion  for  assumed  forms  of  the  time-dependent 
forcing  functions.  There  is  only  a  question  of  the  accuracy  of  the  assumed  form 
for  the  forcing  function  for  the  size  time-step  At  being  used  in  the  analysis.  In 
general,  larger  time-steps  are  likely  to  make  the  assumed  form  for  the  forcing 
function  less  valid. 


1.5  Contents 

Chapter  2  outlines  the  six  algorithms  used  in  this  study  to  compute  the 
response  of  an  SDOF  structural  system.  The  stability  criteria  for  the  three  numer¬ 
ical  integration  methods  and  for  the  numerical  differentiation  method  are  given  in 
Chapter  3.  These  stability  criteria  are  reviewed  and  their  relevance  to  structural 
dynamics/earthquake  engineering  problems  is  discussed.  Also,  a  numerical 
assessment  of  the  largest  time-step  A  tcrWcal  that  can  be  used  in  the  response  anal¬ 
ysis  is  given.  Conclusions  are  made  concerning  the  stability  criteria  for  the  SDOF 
systems  subjected  to  ground  motion  with  At  equal  to  0.02,  0.01  and  0.005  sec¬ 
onds.  A  brief  discussion  and  example  application  of  stability  criteria  for  semidis¬ 
crete  MDOF  structural  system  models  are  also  included. 

Using  damped  SDOF  system  models  with  natural  periods  assigned  based  on 
consideration  of  the  modal  periods  of  hydraulic  structures  providing  significant 
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response  contribution,  an  evaluation  is  made  of  the  accuracy  of  the  computed 
response  values  solved  for  at  regular  time  increments  during  ground  motion.  All 
SDOF  systems  are  assigned  5  percent  damping.  The  results  of  this  extensive 
series  of  numerical  evaluations  are  summarized  in  Chapter  4.  The  results  show 
the  correlation  of  the  accuracy  of  the  six  numerical  step-by-step  procedures,  the 
harmonic  characteristics  of  the  ground  motion,  and  the  time-step  A?  value  (0.02, 
0.01  or  0.005  seconds)  used  in  the  analysis.  Also  included  is  numerical  assess¬ 
ment  of  the  accuracy  of  the  algorithms  for  computing  the  dynamic  response  of 
SDOF  models  in  free  vibration. 

Chapter  5  summarizes  the  results  of  this  study  of  the  accuracy  of  six  numerical 
step-by-step  procedures  used  to  compute  the  dynamic  response  of  SDOF  models 
with  5  percent  damping,  A  brief  discussion  of  the  response  analysis  of  semidis¬ 
crete  MDOF  structural  system  models  to  ground  motion  using  numerical  step-by- 
step  procedures  is  also  included. 

Appendix  A  gives  the  exact  solution  to  a  SDOF  system  subjected  to  a  sine 
wave  base  excitation. 

Appendix  B  gives  the  Fourier  series  representation  for  a  periodic  function  and 
the  response  of  an  SDOF  system  to  a  periodic  force  represented  by  a  Fourier 
series.  This  algorithm  is  in  the  same  category  of  algorithms  as  Duhamel's  integral 
solved  in  a  piecewise  exact  fashion  and  the  Piecewise  Exact  Method  applied 
directly  in  that  it  is  an  exact  solution  to  an  approximation  of  the  actual  loading. 
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2  Six  Numerical  Step-by-Step 
Procedures  of  Analysis  of 
the  Equation  of  Motion  for 
an  SDOF  System 


2.0  Introduction 

This  chapter  outlines  six  numerical  procedures  for  solving  the  dynamic 
response  of  SDOF  models  by  solution  of  the  dynamic  equilibrium  equation  at 
closely  spaced,  discrete  time  intervals  throughout  the  time  of  shaking.  A  base 
acceleration  is  used  for  the  time-dependent  loading.  The  dynamic  response  of  each 
SDOF  system  used  is  characterized  by  the  computed  response  time-histories  of 
accelerations,  velocities,  and  displacements.  The  next  section  begins  with  a  sum¬ 
mary  of  the  equation  of  motion  for  an  SDOF  system  model  subjected  to  a  base 
acceleration  (e.g.,  ground  motion). 

The  six  algorithms  included  in  this  study  are  the  Newmark  |3  method  (with 
values  of  constants  y  and  P  corresponding  to  the  linear  acceleration  method),  the 
Wilson  0  method,  the  Central  Difference  Method,  the  4th  Order  Runge-Kutta 
method,  Duhamel's  integral  solved  in  a  piecewise  exact  fashion,  and  the  Piecewise 
Exact  Method  applied  directly.  All  of  the  algorithms  were  used  in  their  discretized 
forms  (i.e.,  the  loading  and  response  histories  were  divided  into  a  sequence  of  time 
intervals);  thus,  they  are  characterized  as  step-by-step  procedures. 

These  six  algorithms  can  be  categorized  into  two  main  groups,  depending  on 
their  general  approach  to  satisfying  the  differential  equation  of  motion.  The  first 
group  includes  Duhamel's  integral  solved  in  a  piecewise  exact  fashion  and  the 
Piecewise  Exact  Method  applied  directly.  These  two  methods  formulate  exact 
solutions  to  the  equation  of  motion  for  assumed  forms  of  the  time-dependent  forc¬ 
ing  functions  (i.e.,  the  loading  is  approximated  by  a  series  of  straight  lines  between 
the  time-steps).  This  group  is  easily  identified  by  the  fact  that  the  total  response 
consists  of  two  parts:  a  transient  (or  free  vibration)  response  contribution  and  the 
steady-state  response  (or  particular  solution  to  the  specified  form  of  the  loading). 
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The  second  group  of  algorithms  includes  the  Newmark  P  method,  Wilson  0 
method,  Central  Difference  Method,  and  4th  Order  Runge-Kutta  method.  These 
methods  are  referred  to  as  numerical  methods  because  they  approximately  satisfy 
the  equation  of  motion  during  each  time-step  for  the  given  loading.  The  New¬ 
mark  P,  Wilson  0,  and  4th  Order  Runge-Kutta  methods  use  numerical  integration 
to  step  through  the  analysis  of  the  time  response  problem,  where  the  Central  Dif¬ 
ference  Method  uses  numerical  differentiation. 


2.1  Equation  of  Motion  for  SDOF  System 


Consider  the  case  shown  in  Figure  la  of  an  idealized  SDOF  system  subjected 
to  a  time-varying  forcing  function  P(t).  At  time  equal  to  t,  the  SDOF  system  dis¬ 
places  a  distance  x(t)  from  its  at-rest  position  due  to  the  applied  force  of  P(t),  as 
shown  in  Figure  2.  For  a  linear  SDOF  system  acted  on  by  an  externally  applied 
dynamic  force  P(t),  the  equilibrium  criterion  (e.g.,  Chopra  1981  or  Ebeling  1992) 
dictates  that 


fit)  +f(t)  +  fit)  =  Pit)  (1) 

where 

f(t)  =  inertial  force 

f(t)  =  damping  force 

fft)  =  elastic  resisting  force 

P(t)  =  externally  applied  dynamic  force 

The  inertia,  damping,  and  elastic  forces  are  related  to  the  response  quantities  of 
the  SDOF  system  through  the  following  expressions: 


fit)  =  mx(t)  ,  fit)  =  exit)  ,  and  fit)  =  kxit)  (2) 

where 

m  =  mass  of  the  SDOF  system 
xit)  =  relative  acceleration  of  the  mass 
c  =  damping  coefficient  of  the  SDOF  system 
xit)  =  relative  velocity  of  the  mass 
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Figure  1.  Dynamic  response  of  two  damped  SDOF  systems  (Ebeling  1992) 
k  =  stiffness  of  the  SDOF  system 


x{t)  =  relative  displacement  of  the  mass 

Figure  3  shows  that  inertial  force  f  acts  opposite  to  the  acceleration  of  mass  m  at 
time  t. 

Substituting  the  expressions  given  in  Equation  2  into  Equation  1  results  in  the 
equation  of  motion  for  an  SDOF  system: 

mx(t)  +  cx{t)  +  kx{t)  =  P(t)  (3) 

For  earthquake  analyses,  the  dynamic  loading  is  represented  by 

p(  0  =  -Aground*7)  (4) 
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Figure  2.  Forces  acting  on  linear  SDOF  system  at  time  t,  external  force  P(t)  applied  (Ebeling  1992) 


where  5^.ound  (t)  is  the  ground  acceleration  applied  to  the  base  of  the  SDOF  system, 
and  thus  the  equation  of  motion  becomes 

mx(t)  +  cx(t )  +  kx(t )  =  -mXp0UJt)  (5) 

Figure  4  depicts  these  equivalent  dynamic  SDOF  system  problems.  In  alternate 
form,  the  equation  of  motion  may  be  written 

x(t)  +  2  P  cox(  t )  +  co2x(  t)  =  -xvaBBA{  t)  (6) 

where 
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Figure  3.  Inertial  force  acting  opposite  to  the  acceleration  of  mass  m  at  time  t,  external  force  P(t) 
applied  (Ebeling  1992) 


P  =  damping  ratio  =  c/(2mo3) 
(x>  =  circular  frequency  =  <Jklm 
c  =  2mo3fi 


and 


T0  = 


2  TC 


03 


-  2n. 


m 

k 


(7) 


In  earthquake  analyses,  parameters  of  interest  are  relative  displacement,  rela¬ 
tive  velocity,  and  total  acceleration.  The  total  acceleration,  3qotal  (t),  is  simply  the 
sum  of  the  relative  acceleration  plus  the  ground  acceleration 


(8) 
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WO  =  *ground(0  +  *(*) 

=  -{2Po 3x(t)  +  co2x(f)} 
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STEP  1 ,  SOLVE 

mx(t)  +  CX(t)  +  kx(t)  =  -mx  ground  (t) 

STEP  2,  SOLVE 

X  total  W  =  +  Aground  W 


Figure  4.  Equivalent  dynamic  SDOF  system  problems  (Ebeling  1992) 

For  computer  analyses,  the  discretized  form  of  these  equations  is  needed  and 
may  be  represented  by  the  following  notation: 


2.2  Newmark  P  Method 

The  Newmark  P  method  is  based  on  the  following  equations  (e.g.,  Chopra 
1995): 

*,+i  =*,  +  [(1  -  Y)Af]*f  +  (yA t)xj+l  (11) 

and 
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(12) 


*,♦  1  =  *,■  +  (AO*  +  [(0.5  ~P)Af2]x,.  +  (pAf2)x.+1 

where  the  parameters'  p  and  y  define  the  variation  of  response  acceleration  over 
the  time-step  and  control  the  stability  and  accuracy  of  the  method.  Typically,  y  is 
set  equal  to  1/2,  which  corresponds  to  zero  artificial  damping,  and  (3  is  set  to  a 
value  between  1/6  and  1/4.  In  the  analyses  performed  in  this  report  y  -  1/2  and 
P  =  1/6  were  used,  which  corresponds  to  a  linear  variation  of  response  accelera¬ 
tion  over  the  time-step  (i.e.,  the  linear  acceleration  method).  The  original  1959 
version  of  the  Newmark-P  family  of  numerical  methods  required  iteration  to 
implement  Equations  1 1  and  12.  However,  a  modification  can  be  made  to  avoid 
iterations.  This  modified  formulation  is  described  in  this  section. 


Equations  14  and  15  result  from  using  the  following  definitions: 


Ax,  -  x,+i  -x,.  Ax,,  s  x,.+1  -x., 


Ax.  =  x+i  -x.,  and  A Pt  =  PM  ~P, 


(13) 


into  which  Equations  1 1  and  12  are  substituted. 


Ax,.  =  Atx.  +  y  At  Ax', 


Ax,.  =  Atx,.  + 


—  x.  +  PAr2Ax. 
2  ' 


(14) 


(15) 


Solving  Equation  15  for  Ax, 


a  1a  1 

Ax..  =  - Ax..  -  - x. 


1  .. 

■x. 


‘  pAr2  '  PAr  '  2p  ' 
and  then  substituting  into  the  last  term  of  Equation  14  gives 


(16) 


Ax  = 


.  Y 


p  At 


Ax.  - 


-x.  + 
/ 


At 


(17) 


The  incremental  equation  of  motion  (Equation  1 8)  can  be  derived  from  Equa¬ 
tion  3  and  Equation  13: 


m  Ax,  +  c  Ax,.  +  kAxj  =  A P. 


(18) 


1  Note  that  in  Equation  12  and  all  subsequent  equations  in  this  section,  the  variable  (3  describes 
how  the  acceleration  response  of  the  SDOF  system  varies  over  the  time-step  and  does  not  refer  to 
the  damping  ratio,  represented  by  p  in  the  other  sections  of  this  report. 
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Substituting  Equations  16  and  17  into  the  incremental  equation  of  motion  gives 


'  y  1  N 

k  +  c — ■ —  +  m 


pAf  p  At" 


Ax  =  A P  + 


/ 


1  Y 
■m  +  —  c 


l  PA t  P 


2p 


-m  +  At 


(19) 


-X-1 

2p 


Equation  19  can  be  rewritten  as 


k  Ax.  =  AP. 


(20) 


where  k  is  referred  to  as  the  “effective”  stiffness 


k  =  k  +  c- 


+  m- 


1 


P  At  pA  t2 

and  AP  is  referred  to  as  the  “effective”  incremental  force 


(21) 


AP,.  =  AP.  + 


-1—m  +  —  c 

PAt  P  . 


+ 


—  m  +  At 
2P 


'x.,' 


2p 


x.. 


(22) 


Accordingly,  the  incremental  change  in  displacement  Ax,  from  l,  to  h+\  maybe 
determined  by  rearranging  Equation  20  and  from  knowledge  of  the  velocity  and 
acceleration  at  t, 


AP 

-  Ax.  =  — (23) 
k 

Once  Ax,  is  determined,  the  incremental  change  in  velocity  Ax,  and  acceleration 
Ax,  from  t,  to  U+\  may  be  computed  using  Equations  17  and  16,  respectively. 
Rearranging  Equation  13  and  substituting  in  the  values  for  Ax,  and  Ax„  the 
response  velocity  and  acceleration  at  tj+1  can  be  established. 


~  x>;  =  x.  +  Axt  and  xi+I  =  x,  +  Ax.  (24) 

There  are  different  types  of  numerical  methods  in  the  Newmark  P  family 
depending  on  the  values  assigned  to  P  and  y  (Hughes  and  Belytshko  1983;  Hughes 
1987;  Subbaraj  and  Dokainish  1989b;  and  Chopra  1995).  When  the  constant  P  is 
set  equal  to  1/2  and  the  constant  y  is  set  equal  to  1/6,  this  particular  variation  of 
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the  Newmark  (3  family  of  numerical  methods  is  referred  to  as  the  linear  accelera¬ 
tion  method.  The  linear  acceleration  method  is  used  in  the  numerical  studies  to  be 
reported  on  in  subsequent  chapters. 


2.3  Wilson  0  Method 


Although  the  Newmark  p  method  is  versatile  and  accurate,  when  used  with  P 
and  y  values  that  correspond  to  the  linear  acceleration  method,  the  method  is  only 
conditionally  stable,  i.e.,  a  time-step  A t  shorter  than  some  stability  limit  must  be 
used  to  ensure  that  the  solutions  are  bounded1  (e.g.,  Paz  1991  or  Chopra  1995). 
(This  potential  for  instability  is  inherent  in  the  linear  acceleration  method  and  not 
an  artifact  of  the  Newmark  P  method’s  formulation.)  However,  an  uncondition¬ 
ally  stable  form  of  the  linear  acceleration  method  is  the  Wilson  0  method;  for 
0>  1 .37,  the  solution  is  bounded  regardless  of  the  size  of  the  time-step.  The  modi¬ 
fication  that  Wilson  introduced  is  based  on  the  assumption  that  the  response  accel¬ 
eration  varies  linearly  over  an  extended  time  interval  from  t  to  t  +  0  At,  where 
0  >  1 .0.  Note  that  for  0  =  1 .0,  the  Wilson  0  method  is  the  same  as  the  Newmark 
P  method  when  P  =  1/6  and  y  =  1/2. 

Writing  the  equilibrium  criteria  for  an  SDOF  system  at  /,  and  /,  +  6 At,  Equa¬ 
tion  3  becomes 


/,(',)  +/c(',)  (25) 

and 


/,(',  + 0AO  +/c(',-  +  e  At)  +  fk(t.  +  QAt)  =  P(ti  +  Q  At)  (26) 

Subtracting  Equation  25  from  Equation  26  gives 


A f,  +  A fe  +  Afk  =  AP  (27) 

where 

A4  =fI(ti  +  6At)  -/(h)  ,  etc.  (28) 


1  The  stability  limit  for  the  linear  acceleration  method  is  A t  <,  Armucal,  with  A =  0.551(7’J. 
For  the  analyses  of  SDOF  systems  performed  in  this  report,  stability  requirements  are  easily 
satisfied  as  will  be  demonstrated  in  Chapter  3.  However,  for  the  higher  modes  of  vibration  in 
MDOF  systems  or  high-frequency  SDOF  systems,  stability  may  be  an  issue. 
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Accordingly,  the  extended  incremental  form  of  Equation  29  may  be  written  as 


Af.  =  m  Ax  ,  Afc  =  c  Ax  ,  and  Afk  -  k  Ax 


(29) 


where 


Ax;.  =  x(tj  +  QAt)  -  x(tf)  ,  etc. 

and  the  extended  incremental  form  of  the  equation  of  motion  becomes 


(30) 


m  Ax\  +  cAxj  +  kAxi  =  APj  (31) 

Assuming  that  the  response  acceleration  varies  linearly  over  an  extended  time 
interval  from  t  to  t  +  0  A  t. 


x(t)  =  x)  + 


Ax. 

— -  ( t  -  t ,.)  for  t.  <  t  <  tf  +  0Ar 
0Ar  ' 


the  response  velocity  and  displacement  are  given  by: 


(32) 


1  A*,  o 

x(t)  =  x.  +  xdt  -  t)  + - ( t  -  t  )2 

1  2  d At 


and 


x(A  =  x;.  +  xt{t  -  t .)  +  |x,(r  -  tf  +  j^l(t-  f,)3 


(33) 


(34) 


Evaluating  the  response  velocity  (i.e.,  Equation  33)  at  f,  and  +  0A  t  results  in 


x{tt)  =  xt 
and 


(35) 


x{tt  +  0Af)  =  x.  +  x;.0A  /  + 


(36) 
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The  extended  incremental  response  velocity  may  be  obtained  by  subtracting 
Equation  35  from  Equation  36: 


Ax.  =  xflA  t 


+  —  Ax.0At 

2 


(37) 


In  a  similar  fashion,  the  extended  incremental  response  displacement  is  determined 
to  be 


Ax  =  x  0At  + 

l 


-x.(0A t)2  +  -Ax.(0At)2 

2  6 


(38) 


Solving  Equation  38  for  the  extended  incremental  response  acceleration, 


Ax.  = 


(0  At)2 


Ax  - 


(0  At) 


•  x.  -  3x. 


(39) 


and  then  substituting  this  expression  into  the  last  term  of  Equation  37  gives 


a  •  3  i  ,  .  0At  .. 

Ax.  =  - Ax.  -  3x.  -  - x. 

'  0  At  '  '  2  ' 


(40) 


Substituting  Equations  39  and  40  into  the  extended  incremental  equation  of 
motion  (i.e..  Equation  31)  gives 


( 

k  + 

\ 


6 

—  m  + 


0  At 


/ 

+  m 


+  c 


3x.  + 


(41) 


Equation  41  can  be  rewritten  as 


£  Ax.  =  A P 

l  l 


(42) 


where  k  is  again  referred  to  as  the  “effective”  stiffness  but  defined  as 
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j  ,6  3 

k  =  k  +  - m  +  - c 

0  A?  0  A  t 


(43) 


and  AP,  is  again  referred  to  as  the  “effective”  incremental  force  but  defined  as 


aIT  ad  6  ■  o-l  ,•  dAt- 

AP  =  AP.  +  ml  - x.  +  3x  +  cl  3x.  +  - x. 

‘  '  0A t  '  '  '  2  ' 


Accordingly,  the  incremental  change  in  displacement  Ax,  from  t,  to  t,  +  0A?  may 
be  determined  by  rearranging  Equation  42  and  from  knowledge  of  the  velocity  and 
acceleration  at  tv 


Ax.  = 


Once  the  extended  incremental  change  in  displacement  is  known,  Equation  39 
may  be  used  to  compute  the  extended  incremental  change  in  acceleration.  The 
incremental  change  in  acceleration  is  related  to  the  extended  incremental  change  in 
acceleration  through  Equation  46: 


Ax- 

Ax.  =  — - 

'  0 


The  remaining  response  quantities  of  interest  may  be  computed  by  the  following 
expressions: 


Ax.  =  x.At  +  —  Ax  At 
‘  '  2 


Ax.  =  xAt  +  — x.A  t2  +  —  Ax.  At2 
'  '  2  ‘  6  ' 


x,+i  =  +  Ax,. 


*fi  =  *t  + 


=  -(P/+i  -  c^,+i  -  kxi* i) 

wt  '  ' 
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2.4  Centra!  Difference  Method 


The  Central  Difference  Method  is  based  on  the  finite  difference  approximation 
of  the  time  derivatives  of  displacement,  i.e.,  velocity  and  acceleration  (e.g.,  Chopra 
1995).  The  central  difference  expressions  for  velocity  and  acceleration  at  time  t, 
are: 


x. 

I 


'i+1 


xi-i 


2  At 


and 


(52) 


x,.  = 


2  */  +  xi-i 


A  v 


where  at  t,  =  0,  the  initial  response  quantities  x0  and  x0  are  assumed  known1  and  xm] 
is  given  by 


x_,  =  x0  -  A tx0 


(53) 


When  the  expressions  in  Equation  52  are  substituted  in  Equation  3,  the  discretized 
equation  of  motion  may  be  written: 


m 


xi+i 


2x..  +  x 


At1 


M  +  c — — —  +  kx.  =  P 


2  At 


(54) 


or  alternately, 


/ 

m  c 

x.  ,  =  P  - 

( 

m 

\ 

c 

xi-i  ~ 

ft- 2m ) 

KAt2  2  At) 

l+l  l 

K  At2 

2A  t/ 

{  A  t2) 

(55) 


As  with  the  formulations  of  the  Newmark  P  and  Wilson  0  methods,  the  equa¬ 
tion  of  motion  can  be  represented  in  an  analogous  form  to  Hooke’s  law: 


(56) 
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1  The  initial  value  for  relative  acceleration  (x0)  may  be  determined  by  substituting  the  known 
values  of  x0  and  Xj  into  the  equation  of  motion. 
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where  k,  again  the  “effective”  stiffness,  is  written  as 


k  = 


m 

At2 


c 

2  At 


and  P„  again  the  “effective”  force,  is  written  as 


(57) 


c  1 


(58) 


Solving  Equation  56  for  relative  displacement, 


x 


i+i 


r, 

k 


(59) 


Because  the  Central  Difference  Method  is  based  on  the  finite  difference 
approximation  of  the  time  derivatives  of  displacement,  the  determination  of  rela¬ 
tive  velocity  and  acceleration  lags  the  determination  of  relative  displacement  by  At 
(i.e.,  xi+1  is  needed  to  compute  x,  and  x,).  However,  once  xi+1  is  known,  xt  and  x, 
may  be  computed  using  the  expressions  in  Equation  52: 


= 


2  At 


and 


x.  = 


Xi*l  ~2xi+  X,-l 


At7 


2.5  Duhamel's  Integral 


(52  bis) 


Duhamel's  integral  method,  solved  in  a  piecewise  exact  fashion,  idealizes  the 
forcing  function  as  a  succession  of  short-duration  impulses,  with  each  short- 
duration  impulse  being  followed  by  a  free  vibration  response  (e.g.,  Paz  1985; 
Ebeling  1992;  or  Clough  and  Penzien  1993).  Superposition  is  used  to  combine 
each  of  the  short-duration  impulse/ free  vibration  responses  with  the  total  response 
for  the  structural  model.  For  a  continuous  forcing  function,  P{t)  is  divided  into  a 
series  of  pulses  of  duration  dx.  The  change  in  velocity  of  the  SDOF  system  due  to 
the  impulsive  load  may  be  determined  from  Newton’s  law  of  motion: 

m—  —  P(  t)  (60) 

dx 


or 
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(61) 


etc  =  gj > * 
m 

where  P( t)  dx  is  the  impulse,  and  dx  is  the  incremental  velocity.  This  incremental 
velocity  may  be  considered  to  be  an  initial  velocity  of  the  SDOF  system  at  time  x. 

The  solution  to  the  equation  of  motion  for  free  vibration  is 


x(t) 


,  -piof 


X0  +x0Pco 

xQcos(c oDt)  +  -  sin  (c oDt) 


cor 


(62) 


where 


coD  =  G)(/l  -  p2  (63) 

Substituting  Equation  61  for  in  the  second  term  of  Equation  62  and  assuming 
x0  =  0  results  in 

dx(t)  =  e-P“(f--).^(-L^I  sin  <o  (t-T)  (64) 

mcoD 

The  total  relative  displacement  can  be  determined  by  summing  the  differential 
responses,  given  by  Equation  64,  over  the  entire  loading  interval: 


x{t)  =  — - —  rP(x)e'P“('  sin  <x>D(t  -  x)dx  (65) 

mwD  i [ 

Using  the  trigonometric  identity: 


sino>(r-T)  =  sin  cot  cos  cut  -  cos  cot  sin  got 


(66) 


Equation  65  may  be  written 


,  -pcof 


x(t)  =  - \AD{t)  sin  /  -  BJt )  cos  co„rl 

m  co  J 


(67) 


where 


20 


Chapter  2  Six  Numerical  Step-by-Step  Procedures 


(68) 


AD(t )  =  J  P( T)eP“T  cos  go^t  oft 
0 

and 

t 

BD(t)  =f  P( t)cP“t  sin  co0t  dr 
o 

The  expressions  in  Equation  68  can  be  solved  by  several  techniques.  For  this 
study,  the  loading  function  P( r)  is  assumed  to  be  piecewise  linear  and  an  exact 
solution  formulated. 


A P. 

Pi?)  =P(ti. i)  +  —  for  thl  <  t  <  tt 


(69) 


where 


AP,  =  Pit,)  -  P(t ,_,) 


(70) 


When  Equation  69  is  substituted  into  the  expressions  of  Equation  68  and  the  inter¬ 
mediate  variables  Iu  12,  /3,  and  /4  given  in  Equation  71  are  used,  Equation  72 
represents  an  exact  solution. 


,  pO)T 


Ix  =  f  e^x  cos  go0t  o/t  =  - (Poo  cos  oidt  +  o)D  sin  (oDr) 

■'  (Po>)2  +  00  2 


2  r.\  2 

}D 


I2  =  J  e  sin  oonT  dr  = 


eP“- 


-Dw  - (Poo  sin  (oDx  -  u>D  cos  (oDr) 

(poo)2  +coD2 


(71) 


/3  =  f  x  e  P“T  sin  g)dt  dr 


T  - 


Pen 


(pG))2  +G)z)2j 


!  h  + 


G)r 


(Pa)n)2+a)D  9-i 


L  =  J  T  e  COS  G0nT  = 


T  - 


Pg) 


(poo)2  +  C0, 


A  + 


00r 


(pG)^)2  +G)^ 
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where  I{  and  12  are  the  integrals  for  /,  and  I2  before  their  evaluation  at  the  limits. 
By  introducing  the  following  relationships 


AD(tj)  AD{tt_  j)  + 


A P. 


A P, 

J  +  _ L  J 

1  At  4 


*/>(',)  Bd(  t._  j )  + 


A P. 


A P. 
U  +  — 
At 


(72) 


the  relative  displacement,  velocity,  and  acceleration  may  be  determined  using 
Equations  73,  74,  and  75,  respectively. 


X  = 

l 


m  G)r 


AD(tt)  sin  c oDtj  -  BD(t:)  cos 


(73) 


-  p  CO  J1,. 


= 


m<x>r 


{[^(^-pCD^,)]  sin  ( oDtj 

+  [  ^DA d(  0  +  P°)5d(  0  cos  r, } 


(74) 


X,  =  —  (P,  -  cxt  -  kxt) 


(75) 


2.6  Piecewise  Exact  Method 


The  Piecewise  Exact  Method  is  similar  to  the  way  Duhamel’s  integral  was 
solved  in  the  previous  section:  the  forcing  function  is  assumed  to  vary  linearly  in  a 
piecewise  fashion  and  based  on  this  assumption,  an  exact  solution  is  determined 
(Nigam  and  Jennings  1968  and  summarized  in  Appendix  A  of  Gupta  1992). 
However,  the  Piecewise  Exact  Method  is  a  direct  formulation  and  does  not  require 
the  loading  to  be  divided  into  a  series  of  impulses,  as  was  done  in  the  formulation 
of  Duhamel’s  integral. 

Assuming  the  dynamic  loading  varies  in  a  piecewise  linear  fashion: 


A P. 

Pi 0  =  Pit,)  +  -r—it  ~  tj)  for  <  t  <  r  +1 


(76) 
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where 


*P,  =  Pi+ 1  -  Pi 

The  relative  displacement  and  velocity  may  be  determined  by 


(77) 


=  AXi  +  Aground  , 


(78) 


where 


’*/' 

< 

►  x  =  • 

fp-l 

X 

l 

round  i  m 

[p..,J 

(79) 


and 


A  = 


a\\  an 


a2\  a22 


B  = 


bn  bn 
^21  ^22 


(80) 


The  elements  of  the  matrices  A  and  B  are  given  by  the  expressions  in  Equations  8 1 
and  82,  respectively: 


an  =  e 


-poAf 


v  WZ) 


sm  co^A t  +  cos  co^A t 


e 

an  =  -  sin  coDA  t 

G3n 


a)2  e 

a21  - - SU1  “nA? 


0) 


D 


(81) 


_  -pcoAf 
a22  e 


^  A  (OB  •  A  ^ 

cos  o)dA t  -  — —  sm  o)DAt 

C0r-» 


and 
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"n  e 


1  2P2~1  .  p)  sinco^Af 


G)2At  w 


0)r 


/  \ 

2P  +1 

cos  conA  t 

co3  A  t  j 

L/ 

2P 

Cl)3  A  t 


h  =  -e  -P“A t 
"l2  e 


(  \ 
1  2p2  - 1 

lv  °>2A' ; 


sin  coD  A/ 
con 


+  ^P. 


cos  ci)nAr 


co3  A/ 


_  J_  _  2p 
co2  a)3  A/ 


h  -p~  P“A' 
"21  e 


2p2  - 1  |PW 
o)2Ar  w 

/ 

/  2P  ,  1  ^ 


a  o>B  a 

cos  0)D A?  -  — -  sin  (i>DAt 


) 


CO 


V 


(82) 


co3  At  co2 


(coD  sino)DAi  +  Pen  cos  coDAt) 


co2  At 


b 
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-p<oAf 


2  P2  —  1  ' 

co2  A? 


op 


cos  coD  At  -  — -  sin  coD  At 

<*D 


' 
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— —  ( coD  sin  coD  At  +  Po  cos  coDAt) 

Cl)3  A/ 


1 

o2Ar 


Once  the  relative  displacement  and  velocity  have  been  determined,  the  relative 
acceleration  may  be  computed  by 


f/+1  =  ^(P,+  i  "  c*,+i  ~  kxi+i) 


(83) 


2.7  4th  Order  Runge-Kutta  Method 


In  the  application  of  the  4th  Order  Runge-Kutta  method,  the  equation  of  motion 
is  first  reduced  to  two  first-order  differential  equations  (Thomson  1993).  Writing 
the  equation  of  motion  as 


x(t)  =  —  [P(t)  -  kx{t)  -  cx{t)}  =  F(x,x,t)  (84) 

m 
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this  second-order  differential  equation  may  be  written  as  two  first-order  equations: 


x(t)  =  y{t) 
y(t)  =  F(x,y,t) 

Both  x;+1  andyi+1  can  be  expressed  in  terms  of  the  Taylor  series: 


(85) 


*(t;+1)  =  x(tt)  + 

=  y(t,)  + 


^  dx  ^ 


V  dt )  t 


A  t  + 


1  dy 

,dt  St 


A  t  + 


d2x  I  A t1 
V  dt  )t. 


d2y  |  At2 

\  df2 It  2 


(86) 


Ignoring  the  higher  order  derivatives,  and  replacing  the  first  derivative  by  an 
'‘average”  slope,  the  expressions  in  Equation  86  may  be  written: 


x(ti+1)  =  x(tt)  + 

y(h+ 1)  =  yCO  + 


(  dx  ^ 
V  dt  / 

'  dy} 
\  dt ) 


At 


avg 


(87) 


At 


avg 


where,  if  Simpson’s  rule  were  used,  the  “average”  slope  would  be  defined  as 


6 


/ 

+  4 

\ 


dt  /  tt+h/2 


+ 


dy) 


(88) 


In  the  Runge-Kutta  formulation,  the  “average”  slope  is  very  similar  to  that  of 
the  Simpson’s  rule,  except  that  the  center  term  of  Equation  88  is  split  into  two 
terms  and  four  values  of  t,  x,  y,  and  F  are  computed  for  each  point  i  as  follows: 


_£ _ x 

T\  —  tj  X}  =  Xj 


T2  =  t,  +  h/2  X2  =  x,-  +  T,  h/2 

T3  =  t,  +  h/2  X3  =  x,  +  Y2  h/2 

Ta  =  t,  +  h  X4  =  x,  +  Y3  h 


T  =  * 

Y>  =y, 

Y2=yi  +  Flh/2 
Y3  =y,  +  F2  h/2 
Yi,  =y,  +  F3h 


y  _ 

F\  =F(TlM) 


F2  =F(T2x2,Y2) 


f3  =  F(T3X3J3) 
F4  =F(TaXa,Ya) 
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These  quantities  are  then  used  in  the  following  recurrence  equations: 


*  y<V2V2)W,> 

=yM-y,  *  ^(F,+2F2+2F,*Ft) 


(89) 


where  it  is  recognized  that  the  four  values  of  Y  divided  by  6  represent  an  “aver¬ 
age”  slope  dx/dt  and  the  four  values  of  F  divided  by  6  represent  an  “average” 
slope  dy/dt.  Once  the  values  of  the  expressions  in  Equation  89  are  determined,  the 
relative  acceleration  at  time  may  be  computed: 


=  j;{pM  -  -  <*m) 


(90) 
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3  Stability  of  Numerical 
Integration  and  Numerical 
Differentiation  Methods 


3.0  Introduction 


The  stability  criteria  for  the  three  numerical  integration  methods  and  for  the 
numerical  differentiation  method  used  in  the  step-by-step  response  analysis  of  the 
SDOF  systems  analyzed  in  this  study  are  given  in  this  chapter.  Recall  that  the 
numerical  integration  methods  included  in  this  study  are  (a)  the  linear  acceleration 
method  of  the  Newmark  P  family  of  numerical  methods,  (b)  the  Wilson  0  method, 
and  (c)  the  4th  Order  Runge-Kutta  method.  The  numerical  differentiation  method 
used  is  the  Central  Difference  Method.  The  stability  criterion  for  each  of  these 
four  algorithms  is  established  by  the  values  assigned  to  the  constants  that  are  used 
in  the  algorithm  and  the  terms  associated  with  the  structural  model. 

The  stability  condition  requirements  for  numerical  methods  are  categorized  as 
either  unconditional  or  conditional.  A  numerical  method  is  unconditionally  stable 
if  the  numerical  solution  for  any  initial  value  problem  does  not  artificially  grow 
without  bound  for  any  time-step  A  t,  especially  if  the  time-step  is  large  (Bathe  and 
Wilson  1976;  Bathe  1982;  Hughes  and  Belytshko  1983;  Hughes  1987;  and 
Chopra  1995).  The  method  is  conditionally  stable  if  the  previous  statement  is  true 
only  for  those  cases  in  which  the  time-step  A t  is  less  than  some  critical  time-step 
At ^0,1 ■  Figure  5  shows  the  attributes  of  a  stable  response,  computed  using  At  < 
£^t criticab  and  the  attributes  of  a  unstable  response,  computed  using  A t  >  A tcrWcab  for 
the  same  undamped  SDOF  system  in  free  vibration. 
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mx  (t)  +  kx  (t)  =  0 

INITIAL  CONDITIONS:  x(t=0)  =  x0 

X(t=0)  =  x0 


EXACT  SOLUTION:  x  =  x„  cos  CDt  +  sin  CDt 
0  CD 


where:  CD 


-J- 


a. 


t- 

z 

in 

2 

Ui 

O 

5 


a. 

CO 

Q 


time 


Stable  free  vibration  response  with  a  smaller  time-step  than  the  critical 
time-step 


b.  Unstable  free  vibration  response  due  to  a  larger  time-step  than  the 
critical  time-step 


Figure  5.  Example  of  response  for  an  undamped  SDOF  system  in  free 
vibration  (Ebeling  1992) 
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3.1  Stability  Criteria  for  Two  Implicit  Numerical 
Integration  Methods 


The  linear  acceleration  method  of  the  Newmark  P  family  of  numerical  methods 
is  conditionally  stable.  The  critical  time-step  is  defined  as 


At 


critical 


_i _ 1 

n  y Hi  \j  y  -  2  p 


(91) 


and  recall  that  T0  equals  the  natural  (undamped)  period  of  the  SDOF  system 
(Hughes  and  Belytshko  1983;  Hughes  1987;  Subbaraj  and  Dokainish  1989b;  and 
Chopra  1995).  With  y  equal  to  1/2  and  P  equal  to  1/6  for  the  linear  acceleration 
method,  A tcriticai  becomes 


At 


critical 


0.551  T 

O 


Thus,  the  stability  criterion  for  a  SDOF  system  dictates  that 


(92) 


&  *  A?cnhcal 


(93) 


Equation  92  indicates  that  when  the  linear  acceleration  method  is  applied  to  the 
response  analysis  of  SDOF  systems  for  either  free  vibration  or  forced  vibration 
analysis,  the  analysis  requires  two  time-steps  per  natural  vibration  period  of  the 
structure  to  satisfy  stability  criteria.  For  the  case  of  T0  equal  to  0.5  sec,  A tcritical 
becomes  0.276  sec,  and  with  T0  equal  to  0.25  sec,  Apical  reduces  to  0.138  sec. 
Since  all  SDOF  systems  used  in  this  study  are  assigned  T0  equal  to  0.5  sec  or 
0.25  sec  and  are  subjected  to  ground  motion  with  At  set  equal  to  either  0.02,  0.01, 
or  0.005  sec,  it  is  concluded  that  the  numerical  computations  using  the  linear 
acceleration  method  are  stable.  The  results  of  these  forced  vibration  analyses  will 
be  discussed  in  Chapter  4. 

The  Wilson  0  method  is  unconditionally  stable  when  the  value  assigned  to  the 
constant  0  is  greater  than  1.366.  A  value  of  0  equal  to  1 .38  is  used  in  this  study. 
Thus,  no  restraints  (such  as  a  A tcritical  value)  are  placed  on  the  time-step  At  used  in 
the  analyses  from  the  viewpoint  of  numerical  stability  considerations. 
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3.2  Stability  Criteria  for  an  Explicit  Numerical 
Integration  Method 


The  4th  Order  Runge-Kutta  method  is  an  unconditionally  stable  explicit  numer¬ 
ical  integration  method  (Thomson  1993  or  Subbaraj  and  Dokainish  1989b). 
Therefore,  no  restraints  (such  as  a  A tcritical  value)  are  placed  on  time-step  A t  used 
in  the  analyses  from  the  viewpoint  of  numerical  stability  considerations. 


3.3  Stability  Criteria  for  a  Numerical  Differentiation 
Method 

The  Central  Difference  Method  is  a  conditionally  stable  explicit  numerical  dif¬ 
ferentiation  method.  The  critical  time-step  is  defined  as 


aw  ■  -r.  P4) 

Tt 

(Bathe  and  Wilson  1976;  Bathe  1982;  Hughes  and  Belytshko  1983;  Hughes  1987; 
Subbaraj  and  Dokainish  1989b;  Clough  and  Penzien  1993;  and  Chopra  1995). 
Equation  94  indicates  that  when  the  Central  Difference  Method  is  applied  to  the 
response  analysis  of  SDOF  systems  for  either  free  vibration  or  forced  vibration 
analysis,  the  analysis  requires  three  time-steps  per  natural  vibration  period  of  the 
structure  to  satisfy  stability  criteria.  For  the  case  of  Ta  equal  to  0.5  sec,  A tcritical 
becomes  0. 159  sec,  and  with  T0  equal  to  0.25  sec,  A tcritical  reduces  to  0.08  sec. 

Since  the  SDOF  systems  used  in  this  study  (with  T0  equal  to  0.5  sec  or  0.25  sec) 
are  subjected  to  ground  motion  with  A t  set  equal  to  either  0.02,  0.01,  or  0.005  sec, 
it  is  concluded  that  the  numerical  computations  using  the  Central  Difference 
Method  are  stable. 


3.3.1  MDOF  systems 

Stability  conditions  must  be  satisfied  for  each  mode  in  the  MDOF  system 
model,  even  if  the  response  in  the  higher  modes  is  insignificant  (Hughes  1987, 
page  493;  or  Chopra  1995,  page  575).  Accordingly,  stability  criteria,  expressed 
in  terms  of  the  limiting  time-step  A tcriticai  for  conditionally  stable  algorithms,  are 
more  restrictive  for  MDOF  systems  than  for  SDOF  systems.  Hughes  (1987, 
pages  540-542)  shows  a  numerical  exercise  to  establish  the  time-step  A t  value  to 
be  used  in  a  “transient  analysis  of  an  undamped  multidegree-of-ffeedom  structure 
...  for  which  engineering  insight  reveals  that  the  response  will  be  primarily  in  the 
first  six  modes.  Engineering  accuracy  dictates  that  relative  period  error  and 
amplitude  decay  (per  cycle)  be  no  more  than  5  percent  for  any  of  the  first  six 
modes.”  (A  discussion  of  the  issues  related  to  the  accuracy  of  the  numerical  step- 
by-step  procedures  is  postponed  until  Chapter  4.)  However,  this  example  shows 
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that  the  time-step  needed  to  satisfy  the  stability  requirement  of  the  Central  Differ¬ 
ence  Method  is  so  small  that  there  is  no  need  to  worry  about  accuracy  criteria. 

This  is  because  the  maximum  time-step  is  computed  using  Equation  94  with  T0 
replaced  by  the  minimum  period  (i.e.,  the  maximum  frequency)  of  either  the  modes 
for  the  MDOF  model  or  the  individual  elements  modeling  the  structure  (depending 
on  the  formulation  used  to  solve  the  equation  of  motion  of  the  MDOF  system 
model).  The  maximum  time-step  allowed  for  several  unconditionally  stable, 
numerical  step-by-step  procedures  that  may  be  used  for  the  response  analysis  are 
also  computed  in  this  exercise.  The  results  of  Hughes'  example  highlights  the  fact 
that  unconditionally  stable  algorithms  need  be  concerned  only  with  the  issue  of 
accuracy,  and  a  significantly  larger  time-step  A t  may  be  used  in  the  MDOF 
system  response  analysis  being  considered,  compared  to  A tcritical  for  the  Central 
Difference  Method.  Lastly,  Chopra  (1995,  page  170)  observes  that  in  the  analysis 
of  semidiscrete  MDOF  system  models  it  is  often  necessary  to  use  unconditionally 
stable  methods.  Hughes  (1987,  page  536)  notes  that  the  use  of  unconditionally 
stable  methods  is  particularly  important  in  complicated  structural  models 
containing  slender  members  exhibiting  bending  effects. 


3.4  Conclusions 

In  summary,  the  following  conclusions  are  made  regarding  the  stability  require¬ 
ments  of  the  numerical  methods  used  in  this  study: 

a.  The  Wilson  0  method  with  0  equal  to  1 .38  and  the  4th  Order  Runge-Kutta 
method  are  unconditionally  stable  with  no  requirements  made  on  the  time- 
step  At  used  in  the  analyses. 

b.  The  linear  acceleration  method,  of  the  Newmark  P  family,  and  the  Central 
Difference  Method  are  conditionally  stable.  However,  since  the  SDOF 
systems  used  in  this  research  are  subjected  to  ground  motion  with  At  set 
equal  to  either  0.02,  0.01,  or  0.005  sec,  it  is  concluded  that  the  computa¬ 
tions  using  these  two  numerical  methods  are  stable. 

Additionally,  the  following  observation  is  made:  in  most  earthquake 
engineering/dynamic  structural  response  analyses,  a  time-step  At  equal  to  either 
0.02,  0.01,  or  0.005  sec  is  commonly  used  to  define  the  ground  motion  accelera¬ 
tion  time-history.  In  general,  stability  will  not  be  an  issue  for  the  computed  results 
when  either  the  linear  acceleration  method  or  the  Central  Difference  Method  is 
used  for  SDOF  systems  with  T0  ranging  from  0.25  to  0.5  sec.  The  time-step  At 
used  to  accurately  define  the  ground  motion  will  be  much  smaller  than  the  A tcritical 
value  for  either  of  these  numerical  methods.  The  accuracy  of  six  numerical  step- 
by-step  procedures  will  be  discussed  in  detail  in  Chapter  4. 
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4  Accuracy  of  Six  Numerical 
Step-by-Step  Procedures  of 
Analysis  of  the  Equation  of 
Motion  for  SDOF  Systems 


4.0  Introduction 


In  general,  the  accuracy  of  a  numerical  algorithm  is  associated  with  the  rate  of 
convergence  of  the  computed  response  with  the  exact  response  as  At  -  0 
(Hughes  1987).  The  six  algorithms  included  in  this  study  are  the  Newmark  P 
method  (with  values  of  constants  y  and  P  corresponding  to  the  linear  acceleration 
method),  the  Wilson  0  method,  the  Central  Difference  Method,  the  4th  Order 
Runge-Kutta  method,  Duhamel's  integral  solved  in  a  piecewise  exact  fashion,  and 
the  Piecewise  Exact  Method  applied  directly.  Specific  details  regarding  the  equa¬ 
tions  used  in  each  of  the  six  numerical  step-by-step  procedures  are  given  in 
Chapter  2. 

Much  of  the  current  guidance  for  selecting  the  time-step  A  t  used  in  computing 
the  dynamic  response  of  SDOF  and  MDOF  models  to  ground  motion  is  based  on 
studies  of  the  accuracy  of  numerical  methods  for  computing  vibration 
response  of  undamped  SDOF  systems.  The  information  from  the  free  vibration 
studies  is  often  combined  with  useful  but  qualitative  reference  to  the  frequency 
characteristics  of  the  forcing  function.  The  time-step  At  criterion  is  often 
expressed  as  a  fraction  of  the  natural  (undamped)  period  of  the  SDOF  system  for 
a  specified  level  of  accuracy.  Section  4. 1  gives  a  brief  review  of  published  numer¬ 
ical  assessments  of  the  accuracy  of  several  numerical  algorithms  for  different 
time-step  At  values  used  to  compute  the  free  vibration  response  of  undamped 
SDOF  models.  Current  guidance  on  the  factors  to  be  considered  when  choosing 
the  value  of  At  to  be  used  in  response  analysis  of  structures  to  earthquake  shaking 
is  also  included. 

Using  damped  SDOF  system  models  with  natural  periods  assigned  based  on 
consideration  of  the  important  modal  periods  of  hydraulic  structures,  an  evaluation 
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is  made  in  this  study  of  the  accuracy  of  the  computed  response  values  solved  for  at 
regular  time  increments  during  ground  motion.  Recall  from  Section  2. 1  that  a 
ground  acceleration  applied  at  the  base  of  an  SDOF  system  is  equivalent  to  a 
fixed-base  SDOF  system  with  the  forcing  function  applied  directly  to  the  mass. 
Section  4.2  summarizes  the  results  from  an  extensive  series  of  numerical  computa¬ 
tions  used  to  evaluate  the  accuracy  of  the  six  numerical  step-by-step  procedures 
used  in  this  study.  Time-step  A t  values  of  0.02,  0.01,  and  0.005  sec  are  used  in 
the  response  analysis  of  SDOF  systems  subjected  to  base  accelerations  with  differ¬ 
ent  frequency  characteristics.  The  dynamic  response  for  each  damped  SDOF 
structural  model  used  in  this  study  (|3  =  0.05)  is  characterized  by  the  computed 
response  time-histories  of  accelerations,  velocities,  and  displacements.  These 
results,  combined  with  computations  made  using  closed  form  solutions,  allow  for 
the  development  of  quantitative  guidance  as  to  how  the  accuracy  of  the  six  numer¬ 
ical  step-by-step  procedures  is  affected  by  both  the  time-step  At  and  the  fre¬ 
quency  characteristics  of  the  ground  motion. 


4.1  Error  in  Free  Vibration  Response  of  SDOF 
Systems 


The  accuracy  of  a  numerical  step-by-step  procedure  is  usually  characterized 
using  the  computed  results  from  free  vibration  response  analyses  of  undamped 
SDOF  systems  compared  with  the  results  from  a  closed  form  (exact)  solution  to 
the  equation  of  motion.  This  section  briefly  reviews  select  results  from  a  com¬ 
monly  cited  numerical  assessment  of  the  accuracy  of  numerical  algorithms  for 
different  time-step  At  values  used  to  compute  the  dynamic  response  of  SDOF 
structural  models  in  free  vibration.  The  figures  used  to  quantify  the  accuracy  of 
numerical  step-by-step  procedures  and  the  application  example  cited  are  taken 
from  Hughes  (1987). 

Error  is  inherent  in  any  numerical  solution  of  the  equation  of  motion 
(Chopra  1995,  page  170).  A  common  method  used  to  gain  insight  into  the  magni¬ 
tude  of  error  for  a  numerical  step-by-step  procedure  is  to  quantify  the  difference  in 
computed  displacements  with  the  exact  displacements  for  an  undamped  SDOF 
system  in  free  vibration.  The  undamped  SDOF  system  is  set  in  motion  by  an 
initial  displacement  x0  and  an  initial  velocity  x0  at  time  t  =  0  sec.  The  exact 
solution  for  displacement  x  of  the  undamped  SDOF  system  with  time  t  for  this 
boundary  value  problem  is  given  in  Figure  5 .  This  equation  for  x  is  obtained  using 
standard  solution  procedures  for  linear  differential  equations,  such  as  the  method 
of  undetermined  coefficients  (e.g.,  Section  2. 12  in  Kreyszig  1972).  The  response 
described  by  this  equation  and  shown  in  Figure  5  a  is  cyclic  with  a  constant 
maximum  amplitude 


x(t) 


max 


(95) 
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and  constant  circular  frequency  co  ( ( k/m)m  radian).  Recall  from  basic  structural 
dynamics  that  an  undamped  SDOF  system  responds  at  its  natural  (undamped) 
period  T0,  equal  to  2n:/co  sec  (Ebeling  1992).  Thus  the  response  is  said  to  be 
periodic  with  each  complete  response  cycle  occurring  over  a  constant  interval  in 
time  equal  to  T0  sec.  These  attributes  of  a  periodic  response  and  the  same  maxi¬ 
mum  displacement  amplitude  in  each  cycle  facilitate  the  assessment  of  the  error  in 
the  displacements  computed  using  a  numerical  step-by-step  procedure  of  analysis. 

Figure  6  (Hughes  1987)  summarizes  the  results  of  error  assessments  made 
using  several  different  numerical  step-by-step  procedures  to  solve  for  the  dis¬ 
placement  x  of  undamped  SDOF  systems  in  free  vibration  for  a  wide  range  in 
time-step  A  t  values.  The  abscissa  of  Figures  6a  and  6b  is  the  ratio  of  At  divided 
by  Ta,  the  natural  (undamped)  period  of  the  SDOF  system.  Two  definitions  of 
error  are  possible  for  the  free  vibration  problem:  (a)  amplitude  decay,  designated 
as  AD,  and  (b)  period  elongation,  designated  as  PE.  These  two  types  of  errors  are 
shown  in  the  idealized  schematic  in  the  center,  upper  diagram  in  Figure  6  for  one 
complete  cycle  of  displacement  x.  Since  the  SDOF  system  is  undamped,  any 
amplitude  decay  AD  in  displacement  x  (per  cycle)  computed  using  a  numerical 
step-by-step  procedure  will  be  a  measure  of  the  error  in  computed  response.  This 
error  measurement  is  sometimes  reported  as  “algorithmic  damping”  since  the 
actual  response  for  the  SDOF  system  is  undamped  with  no  amplitude  decay  per 
cycle.  Amplitude  decay  is  converted  to  algorithmic  damping  using  the  equation 
given  in  the  center,  upper  schematic.  The  second  type  of  error  possible  is  referred 
to  as  period  elongation  and  measures  the  extension  in  the  time  increment  it  takes  to 
complete  each  cycle  of  harmonic  response. 

The  data  in  Figure  6  show  that  for  the  free  vibration  problem  of  an  SDOF  sys¬ 
tem  with  a  constant  natural  period  Ta,  the  magnitude  of  one  or  both  error  meas¬ 
urements  usually  increases  with  the  time-step  At.  (Refer  to  Hughes  (1987, 

Section  9.3)  for  details  regarding  the  numerical  step-by-step  procedures  identified 
in  this  figure.)  Conversely,  Figure  6  shows  that  for  a  specified  time-step  At,  the 
magnitude  of  one  or  both  error  measurements  is  greater  for  short-period  SDOF 
systems  than  for  long-period  SDOF  systems.  In  summary,  this  figure  shows  the 
errors  associated  with  a  given  numerical  procedure  to  be  a  function  of  (a)  the  time- 
step  At  used  in  the  analysis  and  (b)  the  natural  (undamped)  period  T0.  Similar 
error  plots  are  given  in  Chopra  (1995,  Figure  5.5.2  on  page  173)  and  in  Bathe  and 
Wilson  (1976,  Figure  9.3  on  page  357). 


4.1.1  MDOF  systems 

The  data  given  in  Figure  6  may  also  be  used  to  establish  the  largest  time-step 
At  for  a  specified  level  of  accuracy  in  response  analysis  of  semidiscrete  MDOF 
system  models.  Hughes  (1987,  pages  540-542)  gives  a  numerical  exercise  to 
establish  the  time-step  At  value  to  be  used  in  a  “transient  analysis  of  an  undamped 
multidegree-of-freedom  structure  ...  for  which  engineering  insight  reveals  that  the 
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Figure  6.  Errors  in  free  vibration  response  of  SDOF  systems  for  a-methods,  optimal  collocation  schemes,  and  Houbolt,  Newmark,  Park, 
and  Wilson  methods  (Courtesy  of  Hughes  1987) 


response  will  be  primarily  in  the  first  six  modes.  Engineering  accuracy  dictates 
that  relative  period  error  and  amplitude  decay  (per  cycle)  be  no  more  than  5  per¬ 
cent  for  any  of  the  first  six  modes.”  The  corresponding  value  of  algorithmic 
damping  ratio  is  computed  to  be  0.05/271,  equal  to  0.008.  This  exercise  uses  the 
Figure  6  data  to  establish  the  A t  values  to  be  used  in  each  of  the  numerical  step- 
by-step  procedures  identified  in  this  figure.  From  modal  analysis  of  the  MDOF 
system,  it  is  established  that  the  natural  period  for  the  sixth  mode,  T6,  is  equal  to 
one-tenth  (1/10)  the  natural  period  of  the  first  mode,  T0  (thus,  T0  =  107)). 

Hughes’  exercise  shows  that  when  unconditionally  stable  algorithms,  such  as  the 
Wilson  0  method,  are  used  to  compute  the  response  of  MDOF  system  models,  the 
5  percent  error  criterion  for  amplitude  decay  AD  and  period  elongation  PE 
establishes  two  limiting  values  for  the  ratio  (A t/T6).  For  this  example.  Figure  6 
shows  that  both  limiting  values  of  the  ratio  (A  t/T6)  are  equal  to  0.08.  Note  that 
because  this  is  an  MDOF  system  problem,  the  natural  period  T0  in  the  denomina¬ 
tor  of  the  abscissa  of  Figures  6a  and  6b  is  replaced  by  7),  the  highest  frequency  of 
engineering  significance  contributing  to  system  response  in  this  analysis.  Thus 
the  value  of  the  largest  time-step  A t  that  can  be  used  in  the  response  analysis  using 
the  Wilson  0  method  is  equal  to  0.087)  (0.008 Ta).  Additionally,  if  the  time-step 
criteria  for  AD  and  PE  differ,  the  smallest  value  for  the  ratio  (A t/T6)  is  used  to 
establish  the  largest  time-step  At  since  the  5  percent  maximum  error  criteria  must 
be  satisfied  in  terms  of  both  AD  and  PE. 


4.1.2  Current  guidance  for  assigning  the  time-step  At  to  be  used  in 
earthquake  engineering  dynamic  structural  response  analysis 

Much  of  the  current  guidance  for  selecting  the  time-step  At  used  in  computing 
the  dynamic  response  of  SDOF  and  MDOF  models  to  ground  motion  makes  use  of 
Figure  6  type  data  from  free  vibration  response  analyses  of  undamped  SDOF 
systems.  This  information  is  often  combined  with  useful  but  qualitative  reference 
to  the  frequency  characteristics  of  the  forcing  function.  For  example,  after  going 
through  an  exercise  of  evaluating  the  accuracy  of  numerical  algorithms,  Chopra 
(1995,  pages  172  and  568)  concluded  that  his  Figure  5.5.2  data  (comparable  to 
the  Figure  6  data)  suggest  that  a  time-step  At  equal  to  0.17),,-  would  give  reason¬ 
ably  accurate  results.  (TN  is  the  natural  period  in  seconds  of  the  Nth  mode  of  the 
MDOF  system  model  with  significant  response  contribution.)  This  same  guid¬ 
ance  is  also  given  by  Bathe  and  Wilson  (1976,  pages  351-352),  Clough  and 
Penzien  (1993,  pages  128-129),  and  Paz  (1991,  page  155),  with  the  caveat  of 
when  “the  loading  history  is  simple.”  Chopra  (1995,  pages  172-173)  concludes 
his  discussion  of  computational  error  with  the  observations  that  “...  the  time  step 
should  be  short  enough  to  keep  the  distortion  of  the  excitation  function  to  a  mini¬ 
mum.  A  very  fine  time  step  is  necessary  to  describe  numerically  the  highly  irregu¬ 
lar  earthquake  ground  acceleration  recorded  during  earthquakes;  typically, 

At  =  0.02  seconds  is  chosen  and  this  dictates  a  maximum  time  step  for  computing 
the  response  of  a  structure  to  earthquake  excitation.”  Gupta  (1992,  page  155) 
recommends  that  “the  time  step  used  in  the  response  computations  is  selected  as 
the  smaller  of  the  digitized  interval  of  the  earthquake  acceleragram  or  some  frac¬ 
tion  of  the  period  of  free  vibration,  for  example  7710.” 
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In  summary,  Paz  (1991,  page  155),  Gupta  (1992,  page  155),  Clough  and 
Penzien  (1993,  pages  128-129),  and  Chopra  (1995,  pages  172-173)  all  recognize 
that  the  assignment  of  time-step  A  t  in  response  analysis  to  forced  vibration  should 
consider  the  following  factors:  (a)  the  natural  period  of  the  structure;  (b)  the  rate 
of  variation  of  the  loading  function;  and  (c)  the  complexity  of  the  stiffness  and 
damping  functions.  However,  no  error  summary  similar  to  Figure  6  for  the  forced 
vibration  of  damped  SDOF  systems  in  which  ground  motion  is  represented  by  an 
acceleration  time-history  is  given  in  the  books  on  structural  dynamics  by  Bathe 
and  Wilson  (1976),  Hughes  (1987),  Paz  (1991),  Gupta  (1992),  Clough  and 
Penzien  (1993),  and  Chopra  (1995). 


4.2  Error  in  Response  of  SDOF  Systems  to  Ground 
Motion 


This  section  summarizes  the  results  of  an  assessment  of  the  accuracy  of 
response  of  six  numerical  step-by-step  procedures  used  in  computational  struc¬ 
tural  dynamics.  The  six  algorithms  used  in  this  study  are  representative  of  the 
different  types  of  numerical  procedures  used  to  compute  the  dynamic  structural 
response  to  a  time-dependent  loading.  The  time-dependent  loading  envisioned  in 
this  research  is  that  of  the  motion  of  the  ground  below  a  discrete  structural  model 
and  is  expressed  in  terms  of  a  ground  acceleration  time-history.  The  dynamic 
structural  response  for  each  structural  model  used  in  this  study  is  characterized  by 
the  computed  response  time-histories  of  accelerations,  velocities,  and 
displacements. 

The  six  algorithms  included  in  this  study  are  the  Newmark  P  method  (with  val¬ 
ues  of  constants  y  and  P  corresponding  to  the  lmear  acceleration  method),  the  Wil¬ 
son  0  method,  the  Central  Difference  Method,  the  4th  Order  Runge-Kutta  method, 
Duhamel's  integral  solved  in  a  piecewise  exact  fashion,  and  the  Piecewise  Exact 
Method  applied  directly.  Specific  details  regarding  the  equations  used  in  each  of 
these  numerical  step-by-step  procedures  are  given  in  Chapter  2.  Recall  from  Sec¬ 
tion  2. 1  that  a  ground  acceleration  applied  at  the  base  of  a  SDOF  system  is  equiv¬ 
alent  to  a  fixed-base  SDOF  system  with  the  forcing  function  applied  directly  to  the 
mass. 

These  numerical  results,  combined  with  computations  made  using  closed  form 
solutions,  allow  for  the  development  of  quantitative  guidance  as  to  how  the  accu¬ 
racy  of  the  six  numerical  step-by-step  procedures  are  affected  by  both  the  time- 
step  At  and  the  frequency  characteristics  of  the  ground  motion. 


4.2.1  SDOF  systems 

All  structural  models  used  in  this  numerical  study  are  linear,  SDOF  systems 
with  a  damping  ratio  set  equal  to  5  percent  (P  =  0.05).  Two  SDOF  systems  are 
used  in  this  numerical  study:  T0  =  0.25  sec  (frequency  f0  =  4  Hz)  and  T0  -  0.5  sec 
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{f 0  =  2  Hz).  Recall  from  structural  dynamics  that  frequency/,  (cycles/sec  or  Hz) 
is  equal  to  the  inverse  of  T0  (e.g.,  Ebeling  1992).  The  natural  (undamped)  periods 
and  damping  ratio  of  the  SDOF  systems  used  in  this  numerical  study  were  based 
on  dynamic  response  analyses  procedures  used  to  model  and  analyze  various  types 
of  hydraulic  structures,  such  as  gravity  dams,  arch  dams,  gravity  lock  walls, 
U-ffame  locks,  and  intake  towers.  The  shortest  period  (highest  frequency)  of  engi¬ 
neering  significance  contributing  to  system  response  for  these  hydraulic  structures 
was  also  taken  into  consideration  when  selecting  the  range  in  Ta  values. 


4.2.2  Time-step  A t 

The  time  increments,  At,  used  in  this  numerical  study  are  0.02,  0.01,  and 
0.005  sec.  These  values  are  typical  of  the  At  used  in  discretizing  earthquake 
acceleration  time-histories  recorded  in  the  field  on  strong  motion  accelerographs 
(e.g.,  Hudson  1979). 


4.2.3  Ground  motion 

The  ground  motion  forcing  functions  used  in  this  numerical  study  are  single- 
frequency  harmonics.  The  use  of  a  single  frequency  facilitated  the  evaluation  of 
the  accuracy  of  the  computed  response  values  solved  for  at  regular  time  incre¬ 
ments  during  ground  motion.  Figure  7  shows  the  three  ground  acceleration  time- 
histories  used.  All  three  ground  motions  contain  twenty  cycles  of  sinusoidal 
acceleration  with  peak  ground  acceleration  of  1  g.  The  three  acceleration  time- 
histories  are  distinguished  from  one  another  by  the  time  interval  required  to  com¬ 
plete  each  cycle  of  sinusoidal  acceleration,  designated  as  Tg.  The  values  of  Tg  are 
0.05,  0.25,  and  1.0  sec.  Accordingly,  the  corresponding  durations  of  ground 
motion  are  1,  5,  and  20  sec,  respectively. 

The  three  ground  motions  shown  in  Figure  7,  with  Tg  equal  to  0.05,  0.25,  and 
1  sec,  possess  cyclic  frequencies  fg  of  20  Hz,  4  Hz  and  1  Hz,  respectively.  These 
frequencies  are  often  contained  within  acceleration  time-histories  that  have  been 
recorded  in  the  field  on  strong  motion  accelerographs  during  numerous  earth¬ 
quakes  and  are  often  encountered  in  earthquake  engineering  dynamic  response 
analysis  of  structures.  It  has  also  been  the  experience  of  the  authors  to  encounter 
this  range  of  frequencies  in  the  earthquake  engineering  dynamic  structural 
response  analysis  of  hydraulic  structures. 


4.2.4  Frequency  of  ground  motion  relative  to  frequency  of  SDOF 
systems 

Basic  structural  dynamics  demonstrates  that  the  magnitude  of  the  frequency  fg 
or,  equivalently,  period  Tg,  of  the  forcing  function  relative  to  the  magnitude  of  the 
natural  frequency  f0,  or  period  Ta,  of  the  SDOF  system  impacts  the  magnitude  of 
dynamic  structural  response.  A  response  spectrum  is  a  convenient  plot  for 
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Figure  7.  Equivalent  dynamic  SDOF  system  problems 
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characterizing  these  interrelationships  (Ebeling  1992).  The  response  spectra 
shown  in  Figure  8,  with  j3  =  0.05,  shows  the  relationship  of  the  frequency  content 
for  each  of  the  three  acceleration  time-histories  relative  to  the  natural  frequency  of 
the  two  SDOF  systems.  The  three  ground  motions  used  in  this  study  possess  fre¬ 
quencies  that  are  larger,  equal  to,  and  less  than  the  natural  frequency  of  the  two 
SDOF  systems.  Thus,  three  important  combinations  for  the  frequency  content  of 
the  ground  motion  relative  to  the  natural  frequency  of  the  SDOF  system  are 
included  in  this  investigation.  SDOF  system  responses,  plotted  in  Figure  8  in 
terms  of  pseudo-acceleration  SA  normalized  by  the  peak  ground  acceleration  Amia. 
of  1  g,  demonstrate  that  the  frequency  content  for  the  three  ground  motions  shown 
in  Figure  7  are  sufficiently  close  to  the  natural  frequencies  to  excite  the  SDOF 
systems  and  thus  induce  a  dynamic  response.  (Refer  to  Table  2  in  Ebeling  1992 
for  the  definition  of  SA  and  for  additional  details  regarding  response  spectra.) 


4.2.5  Time-histories  of  432  step-by-step  response  analyses 

The  dynamic  response  for  each  SDOF  model  used  in  this  study  was  character¬ 
ized  by  the  computed  response  time-histories  of  relative  displacements,  relative 
velocities,  relative  accelerations,  and  total  accelerations.  With  three  time-steps 
(A t  =  0.02,  0.01  and  0.005  sec)  three  ground  motions  (Tg  =  0.05,  0.25  and  1  sec), 
a  total  of  432  response  time-histories  were  computed  for  the  two  damped  SDOF 
systems  using  the  six  numerical  step-by-step  procedures.  Additionally,  the  exact 
response  quantities  were  generated  for  each  SDOF  system  for  each  time-step  and 
each  ground  motion  using  the  closed  form  solution  given  in  Appendix  A.  All  cal¬ 
culations  were  made  using  two  computer  programs  developed  for  use  in  this  study. 
Each  of  the  computed  response  time-histories  was  stored  on  disc  in  432  separate 
files.  Each  file  contained  up  to  4,000  response  values  (and  corresponding  time 
t  values),  depending  on  the  time-step  At  value  used  in  the  analysis  and  the  duration 
of  the  ground  motion  used  to  excite  the  damped  SDOF  system. 


4.2.6  Results  from  12  of  the  432  error  studies 

In  this  numerical  study,  the  key  variables  thought  to  impact  the  accuracy  of  the 
results  of  the  six  numerical  step-by-step  procedures  were  (a)  the  time-step  At, 

(b)  the  frequency  content  of  the  ground  motion  (characterized  by  Tg),  (c)  the  value 
of  Tg  relative  to  the  value  of  Ta,  and  (d)  the  natural  period  T0  of  the  SDOF  system. 
The  evaluation  of  the  432  response  time-histories  started  with  a  comparison  with 
their  corresponding  exact  solution.  These  initial  432  comparison  plots  helped  to 
identify  which  variables  contribute  to  the  inaccuracy  in  computed  results  and  the 
extent  of  their  contribution.  Only  12  of  these  432  comparison  figures  are  included 
in  this  report  due  to  space  limitations.  However,  the  time-histories  that  are 
included  in  this  report,  and  discussed  in  the  following  paragraphs,  demonstrate 
how  information  was  extracted  and  used  in  this  extensive  numerical  study. 

Figures  9,  10,  and  1 1  each  show  the  four  response  time-histories  of  an  SDOF 
system  with  T0  equal  to  0.25  sec  and  P  equal  to  0.05.  These  response 
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Figure  8.  Response  spectra  of  two  SDOF  systems  with  5  percent  damping  for  three  harmonic 
forcing  functions 


time-histories  were  computed  using  both  the  Wilson  0  method  (0  =  1.38  and  a 
time-step  A t  equal  to  0.01  sec)  and  the  exact  solution  (Appendix  A).  The  response 
analyses  differ  among  the  three  figures  by  the  frequency  assigned  to  the  ground 
motion,  with  Tg  set  equal  to  0.05  sec  (20  Hz),  0.25  sec  (4  Hz),  and  1  sec  (1  Hz), 
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X  Wilson  e  =  1.38 
-  Exact  solution 


Figure  9.  SDOF  system  [T0  =  0.25  sec)  response  time-histories 

(A t  =  0.01  sec)  computed  using  Wilson  6  =  1.38  for  a  sinusoidal 
forcing  function  of  period  Tg  =  0.05  sec 
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Figure  10.  SDOF  system  ( T0  =  0.25  sec)  response  time-histories 
(Af  =  0.01  sec)  computed  using  Wiison  0  =  1 .38  for  a 
sinusoidal  forcing  function  of  period  T  =  0.25  sec 
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Figure  1 1.  SDOF  system  (T0  =  0.25  sec)  response  time-histories 
(At  =  0.01  sec)  computed  using  Wilson  0  =  1.38  for  a 
sinusoidal  forcing  function  of  period  Tg  =  1 .00  sec 
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respectively.  Visual  comparison  of  the  computed  results  for  each  response  quan¬ 
tity  in  the  three  figures  shows  a  dramatic  variation  in  SDOF  response.  The  rise  in 
amplitude  and  frequency  of  wave  form  for  the  plots  of  relative  displacement  shows 
the  most  dramatic  variation.  Recall  that  these  SDOF  system  responses  are  all 
driven  by  sinusoidal  ground  motions  with  1-g  maximum  amplitude  (Figure  7). 

This  dramatic  variation  in  response  demonstrates  the  impact  that  the  frequency 
content  for  the  ground  motion  relative  to  the  natural  frequency  of  the  SDOF  sys¬ 
tem  has  on  the  dynamic  response.  These  results  also  demonstrate  the  diversity  of 
the  numerical  tests  being  conducted,  even  though  the  forcing  functions  shown  in 
Figure  7  bear  a  strong  resemblance  to  one  another. 

Further  visual  examination  of  Figures  9,  10,  and  1 1  shows  a  larger  error  in 
Figure  9  than  in  the  other  two  figures.  Since  the  frequency  of  the  ground  motion  is 
varied  among  the  three  groups  of  response  analyses  shown  in  Figures  9,  10,  and 
1 1,  the  authors  conclude  that  the  frequency  content  of  the  ground  motion  affects 
the  accuracy  of  all  four  response  parameters  when  using  the  Wilson  0  method.  A 
more  detailed  examination  of  the  results  in  Figure  9  shows  that  the  accuracy  dif¬ 
fers  among  the  four  response  quantities.  For  example,  the  relative  acceleration 
response  data  are  more  accurate  than  the  relative  velocity  response  data. 

Because  these  results  are  all  computed  using  forced  vibration  analyses  of 
damped  SDOF  systems,  the  definitions  of  numerical  errors  used  in  free  vibration 
response  analysis  of  undamped  SDOF  systems  (and  described  in  Section  4.1)  are 
not  appropriate.  A  new  definition  of  numerical  error  in  the  computed  response 
data  is  needed.  Recognizing  that  the  results  from  time-history  analysis  of  linear 
structural  systems  subjected  to  earthquake  excitation  are  usually  concerned  with 
the  extreme  values  in  computed  results,  an  error  definition  is  made  accordingly 
(Figure  12).  The  upper  plot  m  this  figure  is  the  same  relative  displacement  tune- 
history  data  shown  in  Figure  9.  The  time-history  of  100  response  values  (100  = 
duration  of  ground  motion/ A/)  computed  using  the  Wilson  0  method  and  the 
response  values  computed  using  the  exact  method  are  searched  numerically  for 
“peaks  and  valleys.”  This  is  accomplished  using  a  third  computer  program  that 
searches  through  the  response  time-history  data  looking  for  a  reverse  in  sign  in  a 
pair  of  slopes  for  three  adjacent  data  points,  with  the  first  slope  computed  using  a 
pair  of  response  values  at  times  (f,  -  A  t)  and  f,  and  the  next  slope  computed  using 
response  values  at  time  f,  and  (t,  +  At).  The  error  in  the  computed  response  val¬ 
ues,  relative  displacement  in  this  figure,  is  then  computed  for  each  peak  and  each 
valley  in  the  data.  Lastly,  only  those  peaks  and  valleys  with  significant  amplitude, 
say  greater  than  10  percent  of  the  absolute  value  of  the  largest  response  value,  are 
recorded.  An  example  error  calculation  is  made  for  the  first  “peak”  relative  dis¬ 
placement  response  value  and  identified  as  such  in  the  insert  to  the  right  in  Fig¬ 
ure  12.  This  insert  shows  that  the  first  peak  relative  displacement  value  occurs  at 
0.06  sec,  while  the  first  peak  computed  using  the  exact  solution  occurs  at 
0.064  sec.  The  difference  in  time  for  the  two  peaks  attests  to  the  high  frequency  of 
the  response  compared  to  coarse  time-step  used  in  this  Wilson  0  method  response 
analysis  (i.e..  At  =  0.01  sec).  The  computed  error  is  approximately  29  percent. 
This  error  point  and  23  others  are  plotted  versus  time  of  occurrence  (0.06  sec  for 
the  first  peak)  in  the  figure  located  immediately  below  the  relative  displacement 
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forcing  function  of  period  Tg  =  0.05  sec  and  Af  =  0,01  sec 


response  time-history.  This  figure  shows  that  the  error  in  peaks  and  valleys  of 
relative  displacement  ranges  in  value  from  a  low  of  24 . 1  percent  to  a  high  of  3  8 . 4 
percent  throughout  the  duration  of  shaking.  The  maximum  relative  displacement 
for  the  exact  solution  occurs  at  0.042  sec  and  has  a  value  of  -3 .225  mm.  The 
maximum  relative  displacement  computed  by  the  Wilson  0  method  occurs  at 
0.04  sec  and  has  a  value  of  -2.239  mm.  This  corresponds  to  a  30.6  percent  error 
in  maximum  response  computed  by  the  numerical  procedure,  occurring  at  the  first 
“valley”  in  the  response  time-history  (Figure  12). 


4.2.7  Summary  of  numerical  results  from  all  432  error  studies 

An  error  evaluation  similar  to  that  described  in  Section  4.2.6  was  made  for  the 
remaining  431  response  time-histories.  These  error  evaluations  were  performed  on 
all  four  response  variables:  relative  displacement,  relative  velocity,  relative  accel¬ 
eration,  and  total  acceleration.  The  resulting  432  time-history  error  plots  were 
reviewed  by  the  authors.  The  range  in  error  for  all  significant  peaks  and  valleys  of 
a  given  response  parameter  throughout  the  duration  of  shaking  was  tabulated, 
along  with  the  error  in  the  maximum  response  parameter  value.  The  results  of 
these  extensive  error  evaluations  are  summarized  in  Tables  1  through  6  for  relative 
displacement  (designated  Rel.  D),  relative  velocity  (Rel.  V),  relative  acceleration 
(Rel.  A),  and  total  acceleration  (Total  A). 

The  six  algorithms  are  designated  in  Tables  1  through  6  as  follows:  DHM  for 
Duhamel's  integral  solved  in  a  piecewise  exact  fashion;  NMK  for  the  Newmark  P 
method  with  values  of  constants  y  and  (3  corresponding  to  the  linear  acceleration 
method;  PWM  for  the  Piecewise  Exact  Method  applied  directly;  WIL  for  the  Wil¬ 
son  0  method  (0  =  1.38);  CDF  for  the  Central  Difference  Method;  and  RGK  for 
the  4th  Order  Runge-Kutta  method.  Each  table  shares  a  common  value  for  the 
time-step  A t  used  in  the  numerical  analyses.  The  coarsest  time-step  ( At  = 

0.02  sec)  is  used  in  the  numerical  analyses  reported  in  Tables  1  and  2.  The  inter¬ 
mediate  time-step  (0.01  sec)  is  used  in  the  numerical  analyses  reported  in  Tables  3 
and  4.  The  finest  time-step  (0.005  sec)  is  used  in  the  numerical  analyses  reported 
in  Tables  5  and  6.  Each  of  the  three  ground  motions  shown  in  Figure  7  is  distin¬ 
guished  in  the  tables  by  their  Tg  value  (i.e.  0.05,  0.25,  and  1  sec). 

Table  1  summarizes  the  errors  computed  with  At  equal  to  0.02  sec  and  T0 
equal  to  0.25  sec.  The  results  are  presented  in  three  main  groups  and  are  distin¬ 
guished  by  the  ground  motions  used  in  the  analyses.  Recall  that  the  frequency  of 
the  ground  motion  is  reflected  by  the  Tg  value  (and  that  the  ground  motion  fre¬ 
quency^  -  1  /Tg). 

The  results  given  in  Table  1  demonstrate  that  the  accuracy  of  all  six  numerical 
algorithms  depends  on  the  frequency  content  of  the  ground  motion,  with  all  other 
variables  held  constant.  The  magnitudes  of  error  for  all  four  response  parameters 
increase  as  the  frequency  of  the  ground  motion  increases.  Using  a  time-step  At 
equal  to  0.02  sec  for  an  SDOF  system  with  Ta  equal  to  0.25  sec  is  acceptable  for 
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Table  1 

Percentile  Errors  in  Relative  Displacement  (Rel.  D),  Relative  Velocity  (Rel.  V), 

Relative  Acceleration  {Rel.  A),  and  Total  Acceleration  (Total  A)  for  SDOF  System  of 

Tc  =  0.25  sec  and  At  =  0.02  sec 

At 

Error  in  Maximum  Response 
(Range  in  Errors  for  Peak  and  Valley  Values) 

T9,  sec 

Tg 

Parameter 

DHM 

NMK 

PWM 

WIL 

CDF 

RGK 

0.05 

0.4 

Rel.  D 

54 

(50  to  79) 

54 

(52  to  66) 

54 

(50  to  66) 

77 

(70  to  81) 

30 

(40  to  90) 

54 

(53  to  65) 

Rel.  V 

61 

(58  to  87) 

62 

(60  to  85) 

61 

(58  to  86) 

64 

(22  to  80) 

100 

(0  to  1 1 9) 

61 

(58  to  85) 

Rel.  A 

14 

(3  to  45) 

14 

(2  to  44) 

14 

(2  to  45) 

17 

(15  to  46) 

9 

(1  to  41) 

14 

(2  to  43) 

Total  A 

54 

(52  to  66) 

55 

(50  to  67) 

54 

(52  to  66) 

72 

(38  to  89) 

25 

(7  to  91) 

55 

(50  to  66) 

0.25 

0.08 

Rel.  D 

5.3 

(2.2  to  5.3) 

2.2 

(1.3  to  5.1) 

5.3 

(2.2  to  5.3) 

9.9 

(3.8  to  9.9) 

1.9 

(0.1  to  2.9) 

5.4 

(2.4  to  5.5) 

Rel.  V 

3.9 

(2.2  to  4.3) 

3.1 

(1.6  to  4.7) 

3.9 

(2.2  to  4.3) 

11.1 

(4.5  to  11.1) 

4.8 

(1  to  4.8) 

4.4 

(2.2  to  4.4) 

Rel.  A 

5.1 

(1.9  to  5.1) 

0.2 

(0  to  3.1) 

2.9 

(1.9  to  5.1) 

5.6 

(0.3  to  5.6) 

3.1 

(3.1  to  9) 

5.1 

(2  to  5.1) 

Total  A 

4.5 

(2.1  to  4.6) 

2 

(1.3  to  4.2) 

2.6 

(2.1  to  5) 

10.3 

(3.9  to  10.3) 

6.3 

(4  to  7.5) 

4.7 

(2.2  to  5.1) 

1.0 

0.02 

Rel.  D 

0 

0 

0 

0 

0 

0 

Rel.  V 

0 

0 

0 

0 

0 

0 

Rel.  A 

0 

0 

0 

0 

0 

0 

Total  A 

0 

0 

0 

0 

0 

0 

Note:  DHM  =  Duhamel's  integral  solved  in  a  piecewise  exact  fashion. 

NMK  =  Newark  P  method. 

PWM  =  Piecewise  Exact  Method. 

WIL  =  Wilson  0  method. 

CDF  =  Central  Difference  Method. 

RGK  =  4,h  Order  Runge-Kutta  Method. 

long-period  ground  motion  (Tg  =  1  sec,^  =  1  Hz).  However,  for  high-frequency 
ground  motion  (Tg  =  0.05  sec,fg  =  20  Hz)  a  smaller  At  is  required.  The  results 
given  in  Table  2  also  support  this  same  conclusion.  Table  2  differs  from  Table  1 
by  the  value  of  T0  used,  increased  from  0.25  sec  to  0.5  sec. 

Table  3  differs  from  Table  1  in  the  value  of  At  used,  reduced  to  0.01  sec  from 
0.02  sec.  These  results  show  that  the  reduction  in  the  time-step  At  to  0.01  sec 
remarkably  improves  the  accuracy  of  all  six  numerical  methods.  However,  inac¬ 
curacies  are  still  present  in  the  response  values  for  all  six  procedures  for  the 
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Table  2 

Percentile  Errors  in  Relative  Displacement  (Rel.  D),  Relative  Velocity  (Rel.  V), 

Relative  Acceleration  {Rel.  A) 

and  Total  Acceleration  (Total  A)  for  SDOF  System  of 

T0  = 

0.5  sec  and  At  =  0.02  sec 

Error  in  Maximum  Response 

1 

At 

(Range  in  Errors  for  Peak  and  Valley  Values) 

tb 

Parameter 

DHM 

NMK 

PWM 

WIL 

CDF 

RGK 

0.05 

0.4 

57 

57.1 

57 

77.3 

46.5 

57.1 

Rel.  D 

(53.4  to  68.8) 

(55.7  to  72.4) 

(55.5  to  65.4) 

(74.4  to  84.7) 

(19.5  to  89.3) 

(53.5  to  65.3) 

73.9 

74.3 

73.9 

55.2 

74 

73.9 

Rel.  V 

(59.1  to  81.8) 

(59.4  to  93) 

(59.2  to  93.5) 

(2.2  to  104) 

(45.4  to  130) 

(59.2  to  93.5) 

9.7 

9.7 

9.7 

11 

8.3 

9.7 

Rel.  A 

(0.7  to  43) 

(0.7  to  42.9) 

(0.6  to  43) 

(1.1  to  44.4) 

(0.3  to  42.4) 

(0.6  to  43) 

60.8 

60.7 

60.8 

73.4 

45.3 

60.8 

Total  A 

(50.8  to  70.7) 

(51.7  to  74.4) 

(50.8  to  70.7) 

(60.5  to  91.6) 

(16.4  to  77.1) 

(50.8  to  64.8) 

0.25 

0.08 

2.5 

2.8 

2.5 

6.4 

0.1 

2.5 

Rel.  D 

(2.1  to  7.7) 

(0  to  7.0) 

(2.2  to  7.7) 

(4.1  to  25.5) 

(0.1  to  9.2) 

(2.1  to  7.4) 

2.9 

3.5 

2.9 

7.6 

2.2 

2.9 

Rel.  V 

(2.2  to  4.2) 

(1.2  to  5.3) 

(2.2  to  3.7) 

(3.5  to  10.2) 

(0.3  to  4.8) 

(2.2  to  4.2) 

1.7 

2.1 

1.7 

3.8 

2.2 

1.7 

Rel.  A 

(0  to  3.5) 

(0  to  3.8) 

(0  to  3.4) 

(0.2  to  5) 

(0.4  to  4.3) 

(0.2  to  3.5) 

2.2 

2.4 

2.2 

5.9 

4.4 

2.2 

Total  A 

(2  to  5.2) 

(0.6  to  3.7) 

(2.1  to  5.4) 

(1.4  to  18.3) 

(3.4  to  8.3) 

(2  to  5.4) 

1.0 

0.02 

Rel.  D 

0 

0 

0 

0 

0 

0 

Rel.  V 

0 

0 

0 

0 

0 

0 

Rel.  A 

0 

0 

0 

0 

0 

0 

Total  A 

0 

0 

0 

0 

0 

0 

Note: 

DHM 

=  Duhamel's  integral  solved  in  a  piecewise  exact  fashion. 

NMK 

=  Newark  P  method. 

PWM 

=  Piecewise  Exact  Method. 

WIL 

=  Wilson  8  method. 

CDF 

=  Central  Difference  Method. 

RGK 

=  4,h  Order  Runge-Kutta  Method. 

high-frequency  ground  motion  (Tg  =  0.05  sec,fg  =  20  Hz).  The  results  given  in 
Table  4  also  support  this  conclusion.  Table  4  differs  from  Table  3  by  the  value  of 
T0  used,  increased  from  0.25  sec  to  0.5  sec. 

Table  5  differs  from  Table  3  in  the  value  of  A t  used,  reduced  to  0.005  sec  from 
0.01  sec.  These  results  show  that  a  reduction  in  the  time-step  At  to  0.005  sec 
eliminates  all  errors  for  the  six  numerical  methods  in  all  but  the  high-frequency 
ground  motion  compared  with  those  given  in  Table  3.  Small  numerical  inaccu¬ 
racies  are  still  present  in  the  results  for  all  six  numerical  step-by-step  procedures 
for  the  high-frequency  ground  motion  (Tg  =  0.05  sec,fg  =  20  Hz).  The  results 
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Table  3 

Percentile  Errors  in  Relative  Displacement  {Rel.  D),  Relative  Velocity  (Rel.  V), 

Relative  Acceleration  (Rel.  A),  and  Total  Acceleration  (Total  A)  for  SDOF  System  of 

Ta  =  0.25  sec  and  At  =  0.01  sec 

■ 

A  t 

Error  in  Maximum  Response 
(Range  in  Errors  for  Peak  and  Valley  Values) 

m 

Parameter 

DHM 

NMK 

PWM 

WIL 

CDF 

RGK 

0.2 

Rel.  D 

14.8 

(5.3  to  18) 

14.8 

(5  to  20.8) 

13.5 

(6.2  to  19.1) 

30.6 

(24.1  to  38.4) 

7.3 

(7.3  to  21.6) 

13.5 

(12.7  to  18.8) 

Rel.  V 

13.7 

(13.6  to  25.8) 

19.6 

(14.1  to  27.2) 

19.1 

(13.6  to  29.3) 

37.5 

(28.4  to  37.5) 

18.8 

(11.6  to  42.3) 

19.1 

(13.7  to  28.5) 

Rel.  A 

6.6 

(3.1  to  7.5) 

6.65 

(3.1  to  7.5) 

6.6 

(3.1  to  7.5) 

9.3 

(0  to  9.3) 

3.5 

(0.9  to  4.5) 

6.6 

(3  to  7.5) 

Total  A 

14.7 

(8  to  22.5) 

14.9 

(3.1  to  21.1) 

14.7 

(8  to  17.8) 

30.2 

(25.3  to  30.9) 

2.7 

(2.1  to  33.7) 

14.7 

(13.3  to  17.8) 

0.25 

0.04 

Rel.  D 

0 

0 

0 

0 

0 

0 

Rel.  V 

0 

0 

0 

0 

0 

0 

Rel.  A 

0 

0 

0 

0 

0 

0 

Total  A 

0 

0 

0 

0 

0 

0 

1.0 

0.01 

Rel.  D 

0 

0 

0 

0 

0 

0 

Rel.  V 

0 

0 

0 

0 

0 

0 

Rel.  A 

0 

0 

0 

0 

0 

0 

Total  A 

0 

0 

0 

0 

0 

0 

Note:  DHM  =  Duhamel's  integral  solved  in  a  piecewise  exact  fashion. 

NMK  =  Newark  P  method. 

PWM  =  Piecewise  Exact  Method. 

WIL  =  Wilson  0  method. 

CDF  =  Central  Difference  Method. 

RGK  =  4th  Order  Runge-Kutta  Method. 

given  in  Table  6  also  support  this  same  conclusion.  Table  6  differs  from  Table  5 
by  the  value  of  T0  used,  increased  from  0.25  sec  to  0.5  sec. 

The  results  given  in  Tables  1  through  6  show  that  of  the  six  numerical  step-by- 
step  procedures,  the  Wilson  0  method  is  the  least  accurate  and  the  Central  Differ¬ 
ence  Method  is  the  most  accurate.  However,  the  differences  in  the  accuracy  of  the 
computed  results  among  the  six  numerical  step-by-step  procedures  for  the  SDOF 
systems  are  minor  compared  to  the  significance  of  At  and  Tg. 

The  impact  of  the  change  in  natural  period  Ta  from  0.25  sec  to  0.5  sec  for  the 
damped  SDOF  systems  on  the  accuracy  of  the  six  numerical  step-by-step  proce¬ 
dures  is  minor.  Comparison  of  the  results  given  in  Tables  1  and  2  with  Tg  equal  to 
0.25  sec  shows  the  value  of  Tg  relative  to  the  value  of  T0  has  the  most  impact  on 
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Table  4 

Percentile  Errors  in  Relative  Displacement  (Rel.  D),  Relative  Velocity  (Rel.  V), 

Relative  Acceleration  (Rel.  A),  and  Total  Acceleration  (Total  A)  for  SDOF  System  of 

T0  =  0.5  sec  and  At  =  0.01  sec 

At 

Error  in  Maximum  Response 
(Range  in  Errors  for  Peak  and  Valley  Values) 

Tg,  sec 

T9 

Parameter 

DHM 

NMK 

PWM 

WIL 

CDF 

RGK 

0.05 

0.2 

Rel.  D 

14.6 

(12.6  to  16.3) 

14.6 

(10.6  to  16.1) 

14.6 

(12.6  to  15.6) 

30.4 

(28  to  31.9) 

11.9 

(4.8  to  22.4) 

14.6 

(12.6  to  15.6) 

Rel.  V 

13.6 

(13.6  to  27.3) 

13.7 

(13.5  to  52.6) 

13.6 

(13.6  to  37.7) 

30.1 

(29.3  to  50.7) 

13.5 

(13.3  to  25.5) 

13.6 

(13.6  to  52.5) 

Rel.  A 

6 

(3.8  to  6) 

6 

(3.8  to  6) 

6 

(3.8  to  6) 

7.3 

(2.2  to  7.3) 

4.8 

(3  to  4.8) 

6 

(3.8  to  6) 

Total  A 

13.4 

(8.6  to  17.9) 

13.5 

(7.3  to  16.3) 

13.4 

(8.6  to  15.1) 

29.6 

(27.8  to  32.1) 

9.5 

(0.7  to  39.2) 

13.4 

(8.6  to  17.9) 

0.25 

0.04 

Rel.  D 

0 

0 

0 

0 

0 

0 

Rel.  V 

0 

0 

0 

0 

0 

0 

Rel.  A 

0 

0 

0 

0 

0 

0 

Total  A 

0 

0 

0 

0 

0 

0 

1.0 

0.01 

Rel.  D 

0 

0 

0 

0 

0 

0 

Rel.  V 

0 

0 

0 

0 

0 

0 

Rel.  A 

0 

0 

0 

0 

0 

0 

Total  A 

0 

0 

0 

0 

0 

0 

Note:  DHM  =  Duhamel's  integral  solved  in  a  piecewise  exact  fashion. 

NMK  =  Newark  (5  method. 

PWM  =  Piecewise  Exact  Method. 

WIL  =  Wilson  0  method. 

CDF  =  Central  Difference  Method. 

RGK  =  4th  Order  Runge-Kutta  Method. 

accuracy  of  computed  results  when  Tg  and  T0  are  close  to  the  same  value.  How¬ 
ever,  its  influence  is  secondary  compared  to  the  influence  of  A t  and  Tg. 

In  summary,  for  damped  ((3  =  0.05)  SDOF  systems,  the  results  presented  in 
Tables  I  through  6  clearly  demonstrate  that  the  accuracy  of  all  six  numerical 
step-by-step  procedures  depends  primarily  on  (a)  the  value  of  the  time-step  A  t 
used  in  the  response  analysis  and  (b)  the  frequency  content  contained  within  the 
ground  motion.  The  computed  values  of  relative  acceleration  are  slightly  more 
accurate  than  the  computed  values  of  relative  displacement,  relative  velocity,  and 
total  acceleration. 
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Table  5 

Percentile  Errors  in  Relative  Displacement  (Rel.  D),  Relative  Velocity  (Rel.  V), 

Relative  Acceleration  {Rel.  A),  and  Total  Acceleration  (Total  A)  for  SDOF  System  of 

T0  =  0.25  sec  and  At  =  0.005  sec 

n 

At 

(Range 

Error  in  Maximum  Response 
in  Errors  for  Peak  and  Valley  Values) 

Parameter 

DHM 

K9HI 

PWM 

WIL 

CDF 

RGK 

0.1 

Rel.  D 

3.6 

(0.5  to  3.9) 

3.9 

(0.3  to  3.9) 

3.6 

(0.5  to  3.9) 

8.8 

(0.7  to  8.9) 

2 

(0.3  to  6.5) 

3.6 

(0.4  to  3.8) 

Rel.  V 

3.5 

(3.4  to  5.9) 

3.7 

(3.5  to  5) 

3.5 

(3.4  to  5.9) 

9.0 

(8  to  10.6) 

3.9 

(3.2  to  5.9) 

4 

(3.5  to  5.9) 

Rel.  A 

4.7 

(3.5  to  4.8) 

4.7 

(3.7  to  5.1) 

4.7 

(3.6  to  5) 

5.7 

(3.4  to  5.7) 

3.4 

(2.2  to  3.6) 

4.7 

(3.6  to  5) 

Total  A 

4 

(0  to  4) 

4.1 

(0.7  to  4.1) 

4 

(0  to  6.2) 

9.2 

(0  to  9.2) 

3.5 

(0.7  to  11) 

4 

(0  to  8.9) 

0.25 

0.02 

Rel.  D 

0 

0 

0 

0 

0 

0 

Rel.  V 

0 

0 

0 

0 

0 

0 

Rel.  A 

0 

0 

0 

0 

0 

0 

Total  A 

0 

0 

0 

0 

0 

0 

1.0 

0.005 

Rel.  D 

0 

0 

0 

0 

0 

0 

Rel.  V 

0 

0 

0 

0 

0 

0 

Rel.  A 

O 

0 

0 

0 

0 

0 

Total  A 

0 

0 

0 

0 

0 

0 

Note:  DHM  =  Duhamel's  integral  solved  in  a  piecewise  exact  fashion. 

NMK  =  Newark  P  method. 

PWM  =  Piecewise  Exact  Method. 

WIL  =  Wilson  9  method. 

CDF  =  Central  Difference  Method. 

RGK  =  4,h  Order  Runge-Kutta  Method. 

4.2.8  Accuracy  of  numerical  step-by-step  procedures  as  a  function  of 
time-step  Af  and  frequency  contained  within  the  ground  motion 

The  data  contained  in  Tables  1  through  6  show  that  the  accuracy  of  the  six 
numerical  step-by-step  procedures  depends  on  the  value  of  the  time-step  A t  and 
the  frequency  of  the  ground  motion.  This  subsection  describes  two  groups  of  error 
plots  for  each  of  the  four  response  parameters,  given  as  functions  of  the  variables 
At  and  Tg.  Select  data  from  these  tables,  specifically,  the  error  corresponding  to 
the  maximum  responses,  are  reordered  as  a  function  of  the  ratio  of  At  divided  by 
Tg,  in  an  attempt  to  present  this  information  in  a  more  useable  form.  The  first 
group,  Figures  13  through  16,  reports  the  errors  in  relative  displacement,  relative 
velocity,  relative  acceleration,  and  total  acceleration,  respectively,  for  T0  equal  to 
0.25  sec.  The  second  group,  Figures  17  through  20,  reports  the  errors  in  the  four 
response  parameters  for  T0  equal  to  0.5  sec. 
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Table  6 

Percentile  Errors  in  Relative  Displacement  (Rel.  D),  Relative  Velocity  (Rel.  V), 

Relative  Acceleration  (Rel.  A),  and  Total  Acceleration  (Total  A)  for  SDOF  System  of 

Ta  =  0.5  sec  and  At  =  0.005  sec 

A  t 

Error  in  Maximum  Response 
(Range  in  Errors  for  Peak  and  Valley  Values) 

Ts,  sec 

T t 

Parameter 

DHM 

NMK 

PWM 

WIL 

CDF 

RGK 

0.05 

0.1 

Rel.  D 

3.5 

(1.5  to  5) 

3.5 

(0.9  to  5.5) 

3.5 

(2.5  to  5.1) 

8.6 

(7.5  to  9.5) 

2.8 

(2.7  to  5.8) 

3.5 

(3  to  4.3) 

Rel.  V 

3.4 

(3.4  to  4) 

3.4 

(3.4  to  4.1) 

3.4 

(3.4  to  4) 

8.5 

(8.3  to  9.4) 

3.4 

(3.3  to  4.1) 

3.4 

(3.4  to  4.1) 

Rel.  A 

4.6 

(4.3  to  4.9) 

4.4 

(4.3  to  4.9) 

4.4 

(4.3  to  4.9) 

4.9 

(4.1  to  5.1) 

3.9 

(3.8  to  4.4) 

4.4 

(4.3  to  4.9) 

Total  A 

3.3 

(0.8  to  5.6) 

3.3 

(0.7  to  5.2) 

3.3 

(2.5  to  5.6) 

8.4 

(7.5  to  9.1) 

0.1 

(0  to  13) 

3.3 

(2.6  to  5.6) 

0.25 

0.02 

Rel.  D 

0 

0 

0 

0 

0 

0 

Rel.  V 

0 

0 

0 

0 

0 

0 

Rel.  A 

0 

0 

0 

0 

0 

0 

Total  A 

0 

0 

0 

0 

0 

0 

1.0 

0.005 

Rel.  D 

0 

0 

0 

0 

0 

0 

Rel.  V 

0 

0 

0 

0 

0 

0 

Rel.  A 

0 

0 

0 

0 

0 

0 

Total  A 

0 

0 

0 

0 

0 

0 

Note:  DHM  =  Duhamel's  integral  solved  in  a  piecewise  exact  fashion. 

NMK  =  Newark  3  method. 

PWM  =  Piecewise  Exact  Method. 

WIL  =  Wilson  0  method. 

CDF  =  Central  Difference  Method. 

RGK  =  4th  Order  Runge-Kutta  Method. 

These  two  groups  of  data,  Figures  13  through  16  and  Figures  17  through  20, 
demonstrate  that  the  accuracy  of  the  four  response  parameters  computed  using  all 
six  numerical  step-by-step  procedures  clearly  depends  on  the  ratio  of  At  divided 
by  Tg.  Thus,  the  impact  of  the  two  most  important  parameters  on  the  accuracy  of 
the  computed  results  can  be  characterized  in  terms  of  a  single  variable.  Addition¬ 
ally,  the  trends  in  the  data  are  similar  for  each  pair  of  companion  figures,  given  the 
same  response  variable  is  being  considered  (e.g.,  for  Rel.  D  in  Figures  13  and  17; 
Rel.  V  in  Figures  14  and  18;  Rel.  A  in  Figures  15  and  19;  and  Total  A  in  Fig¬ 
ures  16  and  20).  These  eight  figures  show  that  the  error  in  all  response  param¬ 
eters  increases  with  increasing  values  of  the  ratio  A t/Tg.  The  error  in  maximum 
responses  for  the  four  response  parameters  ranges  in  value  from  a  low  of  2.8  per¬ 
cent  to  a  high  of  3 1 .4  percent  with  the  ratio  At/Tg  equal  to  0.2.  However,  reducing 
the  ratio  At/Tg  from  0.2  to  0. 1  reduces  the  error  in  maximum  response  for  all  four 
response  parameters  to  less  than  10  percent. 
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0.02,  0.005  sec 
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Figure  13.  Percentile  errors  in  maximum  relative  displacements  (peak  or  valley  value)  for  SDOF  system  of  T0  =  0.25  sec 


0.005  sec 
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Figure  14.  Percentile  errors  in  maximum  relative  velocity  (peak  or  valley  value)  for  SDOF  system  of  T0  =  0.25  sec 


0.02,  0.005  se 
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Figure  15.  Percentile  errors  in  maximum  relative  acceleration  (peak  or  valley  value)  for  SDOF  system  of  T0  =  0.25  sec 


0.005  sec 


Figure  16.  Percentile  errors  in  maximum  total  acceleration  (peak  or  valley  value)  for  SDOF  system  of  Ta  =  0.25  sec 


0.02,  0.005  sec 
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Figure  17.  Percentile  errors  in  maximum  relative  displacements  (peak  or  valley  value)  for  SDOF  system  of  Ta 


0.02,  0.005  sec 
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Figure  18.  Percentile  errors  in  maximum  relative  velocity  (peak  or  valley  value)  for  SDOF  system  of  T0  =  0.5  sec 


0.005  sec 
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Figure  19.  Percentile  errors  in  maximum  relative  acceleration  (peak  or  valley  value)  for  SDOF  system  of  T0  =  0.5  sec 


0.01  sec 
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Figure  20.  Percentile  errors  in  maximum  total  acceleration  (peak  or  valley  value)  for  SDOF  system  of  Ta  =  0.5  sec 


4.3  Conclusions 


The  accuracy  of  six  numerical  step-by-step  procedures  of  analysis  of  the  equa¬ 
tion  of  motion  for  SDOF  system  models  is  discussed  in  this  chapter.  Current 
guidance  for  the  selection  of  the  time-step  A  t  is  based  on  studies  of  the  accuracy 
of  numerical  methods  using  the  computed  results  from  free  vibration  response 
analyses  of  undamped  SDOF  systems  combined,  frequently,  with  useful  but 
qualitative  reference  to  the  frequency  characteristics  of  the  forcing  function.  This 
time-step  A t  criterion  is  expressed  as  a  fraction  of  the  natural  (undamped)  period 
of  the  SDOF  system  for  a  specified  level  of  accuracy.  The  results  of  the  numeri¬ 
cal  response  calculations  made  in  this  study  of  the  forced  vibration  response  prob¬ 
lem  in  which  a  ground  acceleration  is  applied  at  the  base  of  a  damped  SDOF 
system  adds  to  this  body  of  information. 

Using  damped  SDOF  system  models  with  natural  periods  assigned  based  on 
consideration  of  the  important  modal  periods  of  hydraulic  structures  ( T0  =  0.25 
and  0.5  sec),  an  extensive  numerical  evaluation  is  made  of  the  accuracy  of  the 
computed  response  values  solved  at  regular  increments  in  time  during  ground 
motion.  These  error  assessments  are  given  for  the  four  response  variables  of  rela¬ 
tive  displacement,  relative  velocity,  relative  acceleration,  and  total  acceleration. 
This  assessment  involved  the  evaluation  of  432  response  time-histories.  The 
following  conclusions  are  made  regarding  the  factors  affecting  the  accuracy  of 
results  computed  using  each  of  the  six  numerical  step-by-step  procedures  in  this 
study: 

a.  The  two  variables  shown  to  have  the  most  influence  on  the  accuracy  of  the 
computed  results  for  the  four  response  parameters  are  the  time-step  A t  and 
the  frequency  content  of  the  ground  motion  (characterized  by  f  or, 
equivalently,  7?). 

b.  The  accuracy  in  the  computed  results  for  the  four  response  parameters 
increases  as  the  value  assigned  to  A t  decreases. 

(1)  A  value  of  At  equal  to  0.02  sec  for  the  damped  SDOF  systems  ana¬ 
lyzed  would  be  acceptable  for  long-period  ground  motion  (Tg  =  1  sec, 
fg  =  1  Hz),  but  not  acceptable  for  high-frequency  ground  motion  (Tg 
=  0.05  sec,^  =  20  Hz). 

(2)  A  reduction  in  the  time-step  At  from  0.02  sec  to  0.01  sec  remarkably 
improves  the  accuracy  in  computed  response.  However,  inaccuracies 
are  still  present  in  the  response  values  for  all  six  numerical  step-by- 
step  procedures  for  the  high-frequency  ground  motion  ( Tg  =  0.05  sec, 
fg  =  20  Hz). 

(3)  A  further  reduction  in  the  time-step  At  from  0.01  sec  to  0.005  sec 
eliminates  all  errors  for  the  six  numerical  methods  in  all  but  the  high- 
frequency  ground  motion  (Tg  =  0.05  sec, 7^  =  20  Hz).  The  range  in 
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errors  for  peak  and  valley  values  and  the  error  in  maximum  response 
for  the  four  response  parameters  is  less  than  1 0  percent. 

c.  The  accuracy  in  the  computed  results  for  the  four  response  parameters 
increases  as  the  frequency  content  of  the  ground  motion  decreases  (as  fg  - 
0  or,  equivalently,  as  Tg-  <*>). 

d  The  accuracy  in  the  computed  results  for  the  four  response  parameters  is 
shown  to  correlate  with  the  ratio  A t/Tg.  Thus,  the  impact  of  the  two  most 
important  parameters  on  the  accuracy  of  the  computed  results  can  be  char¬ 
acterized  in  terms  of  a  single  variable.  The  results  show  that  accuracy 
considerations  require  the  value  for  the  ratio  A t/Tg  to  be  less  than  0.2.  A 
value  of  0. 1  for  the  ratio  kt/Tg  is  shown  to  be  sufficiently  accurate  for  all 
six  numerical  step-by-step  procedures. 

e.  The  value  of  Tg  relative  to  the  value  of  T0  has  the  most  impact  on  accuracy 
of  computed  results  when  Tg  and  T0  are  close  to  the  same  value.  However, 
its  influence  is  secondary  compared  with  the  influence  of  A t  and  Tg. 

f.  The  computed  values  of  relative  acceleration  are  slightly  more  accurate 
than  the  computed  values  of  relative  displacement,  relative  velocity,  and 
total  acceleration. 

g.  The  results  show  that  of  the  six  numerical  step-by-step  procedures,  the 
Wilson  0  method  is  the  least  accurate  and  the  Central  Difference  Method 
is  the  most  accurate.  An  improvement  in  the  accuracy  of  results  computed 
using  the  Wilson  0  method  will  be  achieved  if  0  is  set  equal  to  1 .42,  rather 
than  the  1 .38  value  used  in  this  study,  according  to  Chopra  (1995, 

page  581).  However,  the  differences  in  the  accuracy  of  the  computed 
results  among  the  six  numerical  step-by-step  procedures  for  the  SDOF  sys¬ 
tems  are  minor  compared  with  the  significance  of  At  and  Tg. 
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5  Results  and  Conclusions 


5.0  Introduction 


This  report  summarizes  an  assessment  of  the  accuracy  of  six  numerical  step- 
by-step  procedures  used  in  computational  structural  dynamics.  The  six  algorithms 
used  in  this  study  are  representative  of  the  different  types  of  procedures  used  to 
compute  the  dynamic  structural  response  to  a  time-dependent  loading.  The  time- 
dependent  loading  envisioned  in  this  research  is  that  of  the  motion  of  the  ground 
below  a  discrete  structural  model  and  is  expressed  in  terms  of  a  ground  accelera¬ 
tion  time-history.  The  dynamic  structural  response  for  each  structural  model  used 
in  this  study  is  characterized  by  the  computed  response  time-histories  of  accelera¬ 
tions,  velocities,  and  displacements. 

All  structural  models  used  in  this  study  were  linear,  single-degree-of-ffeedom 
(SDOF)  systems.  The  natural  (undamped)  periods  T0  of  these  SDOF  systems 
were  selected  based  on  consideration  of  the  important  modal  periods  of  hydraulic 
structures  such  as  gravity  dams,  arch  dams,  gravity  lock  walls,  U-ffame  locks,  and 
intake  towers.  Each  of  the  forcing  functions  used  in  this  study  was  single¬ 
frequency  harmonics.  The  use  of  a  single  frequency  facilitated  the  evaluation  of 
the  accuracy  of  the  computed  response  values  solved  for  at  regular  time 
increments  during  ground  motion. 

The  time  increments  A t  used  in  this  study  were  0.02,  0.01,  and  0.005  sec. 

These  values  are  typical  of  the  A t  used  in  discretizing  earthquake  acceleration 
time-histories  recorded  in  the  field  on  strong  motion  accelerographs. 

The  six  algorithms  included  in  this  study  were  the  Newmark  P  method  (with 
values  of  constants  y  and  P  corresponding  to  the  linear  acceleration  method),  the 
Wilson  0  method,  the  Central  Difference  Method,  the  4th  Order  Runge-Kutta 
method,  Duhamel's  integral  solved  in  a  piecewise  exact  fashion,  and  the  Piecewise 
Exact  Method  applied  directly.  All  of  these  algorithms  were  used  in  their  discre¬ 
tized  forms  (i.e.,  the  loading  and  response  histories  were  divided  into  a  sequence  of 
time  intervals);  thus,  they  are  characterized  as  step-by-step  procedures. 

The  selection  of  the  size  of  the  time-step  At  to  be  used  in  the  step-by-step  cal¬ 
culation  of  the  dynamic  response  of  the  SDOF  (and  of  MDOF  semidiscrete 
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structural  models)  is  restricted  by  stability  and/or  accuracy  considerations.  The 
primary  requirement  of  a  numerical  algorithm  is  that  the  computed  response  con¬ 
verge  to  the  exact  response  as  At  -  0  (Hughes  1987).  However,  the  number  of 
computations  increases  as  the  time-step  A t  is  made  smaller  in  a  dynamic  analysis, 
an  important  issue  for  response  analysis  of  semidiscrete  MDOF  structural  system 
models. 

In  addition  to  accuracy  considerations,  stability  requirements  must  also  be  con¬ 
sidered  during  the  selection  of  the  time-step  At  to  be  used  in  a  step-by-step 
response  analysis  by  either  of  the  three  numerical  integration  methods  or  by  the 
numerical  differentiation  method.  Stability  criterion  is  expressed  in  terms  of  a 
maximum  allowable  size  for  the  time-step  A tcritical.  The  value  for  AtcnticaI  differs 
among  the  four  numerical  algorithms. 

No  stability  criterion  (expressed  in  terms  of  a  limiting  time-step  value)  is 
needed  for  Duhamel's  integral  solved  in  a  piecewise  exact  fashion  and  the  Piece- 
wise  Exact  Method  applied  directly.  These  two  methods  formulate  exact  solutions 
to  the  equation  of  motion  for  assumed  forms  of  the  time-dependent  forcing  func¬ 
tions.  There  is  only  a  question  of  the  accuracy  of  the  assumed  form  for  the  forcing 
function  for  the  size  time-step  At  being  used  in  the  analysis.  In  general,  larger 
time-steps  are  likely  to  make  the  assumed  form  for  the  forcing  function  less  valid. 


5.1  Stability  Requirements  for  the  Four  Numerical 
Methods  Used  for  Response  Analysis 


Stability  requirements  for  numerical  methods  are  expressed  in  terms  of  a  lim¬ 
iting  or  critical  time-step  A tcritical.  The  stability  criteria  for  the  three  numerical 
integration  methods  and  for  the  numerical  differentiation  method  used  in  this  study 
are  evaluated  in  Chapter  3.  The  following  conclusions  are  made  regarding  the 
stability  requirements: 

a.  The  Wilson  0  method  with  0  equal  to  1 .38  and  the  4th  Order  Runge-Kutta 
method  are  unconditionally  stable  with  no  requirements  made  on  the  time- 
step  At  used  in  the  analyses. 

b.  The  linear  acceleration  method,  of  the  Newmark-P  family,  and  the  Central 
Difference  Method  are  conditionally  stable.  However,  since  the  SDOF 
systems  used  in  this  research  are  subjected  to  ground  motion  with  At  set 
equal  to  either  0.02,  0.01  or  0.005  sec,  it  is  concluded  that  the  computa¬ 
tions  using  these  two  numerical  methods  are  stable. 

Additionally,  the  following  observation  is  made:  in  most  earthquake 
engineering/dynamic  structural  response  analyses,  a  time-step  At  equal  to  either 
0.02,  0.01,  or  0.005  sec  is  commonly  used  to  define  the  ground  motion  accel¬ 
eration  time-history.  In  general,  stability  will  not  be  an  issue  for  the  computed 
results  when  using  either  the  linear  acceleration  method  or  the  Central  Difference 
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Method  for  SDOF  systems  with  T0  ranging  from  0.25  sec  to  0.5  sec.  The  time- 
step  At  used  to  accurately  define  the  ground  motion  will  be  much  smaller  than  the 
At criticai  value  for  either  of  these  numerical  methods. 


5.2  Accuracy  of  Response  of  SDOF  Systems  to 
Ground  Motion 


A  review  is  made  in  Section  4.1.2  of  the  current  guidance  for  selecting  the 
time-step  At  used  in  computing  the  dynamic  response  of  SDOF  and  MDOF 
models  to  ground  motion.  Recall  that  a  ground  acceleration  applied  at  the  base  of 
an  SDOF  system  is  equivalent  to  a  fixed-base  SDOF  system  with  the  forcing  func¬ 
tion  applied  to  the  mass.  Current  guidance  is  shown  to  be  based  on  studies  of  the 
accuracy  of  numerical  methods  using  the  computed  results  from  free  vibration 
response  analyses  of  undamped  SDOF  systems  combined,  frequently,  with  useful 
but  qualitative  reference  to  the  frequency  characteristics  of  the  forcing  function. 
This  criterion  is  often  expressed  as  a  fraction  of  the  natural  (undamped)  period  of 
the  SDOF  system  for  a  specified  level  of  accuracy.  The  current  guidance  given  by 
Chopra  (1995),  Gupta  (1992),  and  the  American  Society  of  Civil  Engineers 
(1986)  Standard  ASCE  4-86,  on  the  factors  to  be  considered  in  choosing  the  value 
of  At  is  as  follows: 

a.  Chopra  (1995,  pages  172-173)  concludes  his  discussion  of  computational 
error  with  the  observations  that  “...  the  time  step  should  be  short  enough  to 
keep  the  distortion  of  the  excitation  function  to  a  minimum.  A  very  fine 
time  step  is  necessary  to  describe  numerically  the  highly  irregular  earth¬ 
quake  ground  acceleration  recorded  during  earthquakes;  typically,  At  = 

0.02  seconds  is  chosen  and  this  dictates  a  maximum  time  step  for  comput¬ 
ing  the  response  of  a  structure  to  earthquake  excitation.” 

b.  Gupta  (1992,  page  155)  recommends  that  “the  time  step  used  in  the 
response  computations  is  selected  as  the  smaller  of  the  digitized  interval  of 
the  earthquake  acceleragram  or  some  fraction  of  the  period  of  free  vibra¬ 
tion,  for  example  T/10.” 

c.  The  ASCE  4-86  Standard  (page  21)  states  that  an  acceptable  rule  for  time- 
history  response  analysis  is  that  the  time-step  At  used  be  small  enough 
such  that  the  use  of  1/2A t  does  not  change  the  response  by  more  than 

10  percent.  For  the  commonly  used  numerical  step-by-step  procedures  of 
Wilson  0  and  Newmark  P,  the  maximum  time-step  size  should  be  one- 
tenth  (1/10)  the  shortest  period  of  interest,  and  for  the  Piecewise  Exact 
Method  (or  Nigam-Jennings  Method)  the  maximum  time-step  size  should 
be  one-fifth  (1/5)  the  shortest  period  of  interest.  Normally,  the  shortest 
period  of  interest  need  not  be  less  than  0.03  sec  (33  Hz  for  nuclear 
structures). 
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Pa lz  (1991,  page  155),  Gupta  (1992,  page  155),  Clough  and  Penzien  (1993, 
pages  128-129),  and  Chopra  (1995,  pages  172-173)  all  recognize  that  the  assign¬ 
ment  of  time-step  A  t  in  response  analysis  to  forced  vibration  should  consider  the 
following  factors: 

a.  The  natural  period  of  the  structure. 

b.  The  rate  of  variation  of  the  loading  function. 

c.  The  complexity  of  the  stiffness  and  damping  functions. 

Using  damped  SDOF  system  models  ((3  =  0.05)  with  natural  periods  assigned 
based  on  consideration  of  the  modal  periods  of  hydraulic  structures  with  signifi¬ 
cant  response  contribution  ( T0  =  0.25  and  0.5  sec),  an  extensive  numerical  evalua¬ 
tion  is  made  of  the  accuracy  of  the  computed  response  values  solved  for  at  regular 
increments  in  time  during  ground  motion.  The  numerical  results,  given  in  Tables  1 
through  6  and  in  Figures  13  through  20,  show  the  interrelationship  between  the 
accuracy  of  the  six  numerical  step-by-step  procedures  with  the  harmonic  charac¬ 
teristics  of  the  ground  motion  (using  Tg  =  0.05,  0.25  and  1  sec)  and  the  time-step 
At  value  (0.02,  0.01  or  0.005  sec)  used  in  the  analysis  of  each  of  the  damped 
SDOF  systems.  These  error  assessments  are  given  for  the  four  response  variables 
of  relative  displacement,  relative  velocity,  relative  acceleration,  and  total  accel¬ 
eration.  The  following  conclusions  are  made  regarding  the  accuracy  of  the  six 
numerical  step-by-step  procedures  used  in  this  study: 

a.  All  six  numerical  step-by-step  procedures  provide  accurate  results  for 
ground  motion  with  a  range  in  frequency  fg  from  1  Hz  to  20  Hz  (or,  equiv¬ 
alently,  Tg  from  0.05  sec  to  1  sec)  when  a  0.005-sec  time-step  is  used  in 
the  numerical  analysis. 

b.  The  accuracy  of  computed  results  diminishes  with  increasing  time-step 
value  used  in  the  numerical  step-by-step  analysis. 

c.  Numerical  errors  are  observed  when  solving  for  the  SDOF  system 
response  to  high-frequency  ground  motion  (fg  =  20  Hz  or  Tg  =  0.05  sec) 
when  the  time-step  At  is  increased  from  0.005  sec  to  0.01  sec  in  the 
numerical  response  analysis.  No  numerical  errors  are  observed  when 
solving  for  the  SDOF  system  response  to  intermediate  and  low-frequency 
ground  motions  (fg-  4  and  1  Hz  or  Tg  =  0.25  and  1  sec)  when  the  time- 
step  At  is  set  equal  to  0.01  sec  in  the  numerical  analysis. 

d.  Numerical  errors  are  observed  when  solving  for  the  SDOF  system 
response  to  the  entire  spectrum  of  ground  motion  frequencies  (fg  =  20  Hz 
or  Tg  =  0.05  sec ,fg  =  4  Hz  or  Tg  =  0.25  sec,  and  fg  =  1  Hz  or  Tg  -  1  sec) 
when  the  time-step  At  is  increased  from  0.01  sec  to  0.02  sec  in  the  numer¬ 
ical  analysis.  The  magnitudes  of  these  errors  increase  as  the  frequency  fg 
contained  within  the  ground  motion  increases.  The  magnitude  of  error  for 
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some  response  parameters  is  likely  to  be  unacceptable  with  At  set  equal  to 
0.02  sec  in  the  response  analysis. 

e.  The  accuracy  in  the  computed  results  for  the  four  response  parameters 
increases  as  the  frequency  content  of  the  ground  motion  decreases  (as  fg  - 
0  or,  equivalently,  as  Tg  -  °°). 

f.  The  accuracy  in  the  computed  results  for  the  four  response  parameters  is 
shown  to  correlate  with  the  ratio  A t/Tg.  Thus,  the  impact  of  the  two  most 
important  parameters  on  the  accuracy  of  the  computed  results  can  be  char¬ 
acterized  in  terms  of  a  single  variable.  The  results  show  that  accuracy 
considerations  require  the  value  for  the  ratio  At/Tg  to  be  less  than  0.2.  A 
value  of  0. 1  for  the  ratio  At/Tg  is  shown  to  be  sufficiently  accurate  for  all 
six  numerical  step-by-step  procedures. 

g.  The  value  of  Tg  relative  to  the  value  of  T0  has  the  most  impact  on  accuracy 
of  computed  results  when  Tg  and  T0  are  close  to  the  same  value.  However, 
its  influence  is  secondary  compared  with  the  influence  of  At  and  Tg. 

h.  The  computed  values  of  relative  acceleration  are  slightly  more  accurate 
than  the  computed  values  of  relative  displacement,  relative  velocity,  and 
total  acceleration. 

i.  The  results  show  that  of  the  six  numerical  step-by-step  procedures,  the 
Wilson  0  method  is  the  least  accurate  and  the  Central  Difference  Method 
is  the  most  accurate.  The  accuracy  of  results  computed  using  the  Wilson 
0  method  will  be  improved  if  0  is  set  equal  to  1 .42,  rather  than  1.38  used 
in  this  study,  according  to  Chopra  (1995,  page  581).  However,  the  differ¬ 
ences  in  the  accuracy  of  the  computed  results  among  the  six  numerical 
step-by-step  procedures  for  the  SDOF  systems  are  minor  compared  with 
the  significance  of  At  and  Tg. 

The  authors  envision  that  the  error  summaries  given  in  Tables  1  through  6  and 
in  Figures  13  through  20  may  be  used  to  establish  an  acceptable  time-step  At 
value  to  be  used  in  earthquake  engineering  dynamic  response  analyses  of  SDOF 
structures  (with  P  =  0.05)  using  any  one  of  the  six  numerical  step-by-step  pro¬ 
cedures.  The  scenario  may  be  as  follows:  Information  regarding  the  frequencies 
contained  within  the  acceleration  time-history  is  made  available  (e.g.,  from 
response  spectra  or  Fourier  amplitude  plots).  For  a  predefined  level  of  accuracy 
in  response  parameter(s)  and  given  knowledge  of  the  frequency  characteristics  of 
the  ground  motion  (i.e.,fg  or  Tg),  these  tables  and  figures  may  be  used  to  establish 
the  minimum  At  value  for  accuracy  considerations.  The  three  possible  outcomes 
are  as  follows: 

a.  The  minimum  At  value  for  accuracy  considerations  is  smaller  than  the  At 
value  contained  within  this  ground  motion  that  was  initially  proposed  for 
use  in  the  structural  response  analysis.  This  situation  will  require  replace¬ 
ment  of  the  initial  ground  acceleration  time-history  (used  to  generate  this 
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initial  Tg  infonnation)  with  one  that  has  been  processed  at  this  finer 
time-step  so  that  numerical  error(s)  contained  within  the  dynamic  struc¬ 
tural  response  analysis  will  be  within  acceptable  levels. 

b.  The  minimum  At  value  for  accuracy  considerations  is  much  larger  than  the 
At  value  contained  in  this  (initial)  acceleration  time-history.  Considera¬ 
tions  related  to  computational  efficiency  in  the  dynamic  structural  response 
analysis  (i.e.,  engineering  costs)  may  dictate  that  a  coarser  time-step  be 
used  in  the  response  analysis.  The  time-step  may  be  increased  by  first 
transferring  the  acceleration  time-history  signal  to  the  frequency  domain 
and  then  returning  the  signal  back  to  the  time  domain  but  with  a  new, 
coarser  time-step. 

c.  The  minimum  At  value  for  accuracy  considerations  is  equal  to  the  At  value 
contained  in  the  initially  selected  ground  motion.  The  dynamic  structural 
response  analysis  then  proceeds  using  this  ground  motion. 

5.3  Baseline  Correction  of  Ground  Motion 

As  a  general  rule  to  avoid  inaccurate  response  predictions,  the  ground  accel¬ 
eration  time-history  used  to  represent  earthquake  excitation  in  the  dynamic 
response  analyses  of  linear  SDOF  and  semidiscrete  MDOF  structural  models  shall 
be  baseline  corrected.  Baseline  correction  is  essential  in  the  response  analysis  of 
nonlinear  structural  systems. 


5.4  Observations  Made  Regarding  Response  Anal¬ 
ysis  of  Semidiscrete  MDOF  System  Models 


Discussions  in  this  report  regarding  the  accuracy  of  computed  response  using 
six  numerical  step-by-step  procedures  have  focused  on  the  response  analysis  of 
SDOF  systems.  This  section  discusses  some  observations  made  regarding 
response  analysis  of  semidiscrete  MDOF  system  models. 

In  general,  the  numerical  step-by-step  procedures  used  in  solving  linear  SDOF 
systems  in  this  report  are  easily  extended  to  deal  with  MDOF  models  by  replacing 
the  Chapter  2  scalar  quantities  by  matrices.  This  type  of  formulation  is  a  direct 
solution  of  the  equation  of  motion  (the  matrix  form  of  Equation  5  in  Chapter  2). 
Thus,  the  solution  for  the  time-history  of  response  is  performed  directly  by  a 
numerical  step-by-step  algorithm  dealing  simultaneously  with  all  degrees  of  free¬ 
dom  (DOF)  in  the  response  vector.  However,  for  an  MDOF  model  of  a  hydraulic 
structure,  such  as  an  arch  dam,  with  a  large  number  of  DOF's,  it  is  computa¬ 
tionally  advantageous  to  transform  the  equation  of  motion  to  modal  coordinates 
before  carrying  out  the  time  response  analysis.  The  reason  is  that,  in  most  cases, 
the  significant  response  of  the  dam  structure  can  be  adequately  described  by  the 
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few  lowest  vibration  modes,  and  thus  solution  of  the  complete  set  of  equations  is 
avoided  (Ghanaat  1993).  For  example.  Engineer  Manual  (EM)  1 1 10-2-2201 
(Headquarters,  U.S.  Army  Corps  of  Engineers,  1994,  pages  7-12  and  7-13)  states 
that  for  linear-elastic  response,  a  sufficient  number  of  modes  should  be  included  so 
that  at  least  90  percent  of  the  total  dynamic  response  is  achieved.  For  large  arch 
dams  this  usually  involves  all  vibration  modes  with  frequencies  less  than  10  Hz. 
Smaller  arch  dams,  which  are  stiffer,  may  reach  their  90  percent  level  of  responses 
when  all  vibration  modes  under  20  Hz  are  considered. 

This  second  type  of  formulation  makes  use  of  modal  methods  to  transform  the 
extensive  number  of  equations  of  motion  for  the  MDOF  system  model  into  uncou¬ 
pled  modal  equations  for  a  much  smaller  series  of  SDOF  equations  of  motion. 
These  transformed  SDOF  equations  of  motion  are  solved  at  each  step  in  time 
(i.e.,  time  t,  t  +  At,  t  +  2A t,  etc.).  Superposition  is  used  to  combine  responses 
computed  by  each  SDOF  equation  of  motion  in  each  mode,  for  the  complete 
dynamic  response  of  MDOF  model  at  each  increment  in  time.  A  change  from  the 
actual  (i.e.,  finite  element)  coordinate  basis  is  made  to  the  basis  of  eigenvectors  for 
the  modal  generalized  equations  in  this  formulation.  Lastly,  because  superposition 
is  employed  in  the  analytical  formulation,  use  of  this  procedure  is  restricted  to 
linear  MDOF  models.  Refer  to  Bathe  and  Wilson  (1976),  Clough  and  Penzien 
(1993),  Ghanaat  (1993),  or  Chopra  (1995)  for  additional  details  regarding  this 
formulation. 


5.4.1  Stability  requirements 

The  stability  criterion  for  the  linear  acceleration  method  and  the  Central  Differ¬ 
ence  Method,  expressed  in  terms  of  a  critical  time-step  A tcritical,  was  shown  not  to 
be  restrictive  on  the  response  analyses  for  the  SDOF  systems  analyzed  in  Chap¬ 
ter  4.  Recall  that  for  the  SDOF  system  with  T0  equal  to  0.25  sec,  the  smallest 
values  of  A tcriticai  are  computed  to  be  0. 138  sec  for  the  linear  acceleration  method 
and  0.08  sec  for  the  Central  Difference  Method.  The  time-steps  A t  used  in  the 
SDOF  response  analyses  are  set  equal  to  0.02,  0.01,  and  0.005  sec  (Chapter  4). 

The  stability  of  a  numerical  method  is  a  critical  consideration  in  the  analysis  of 
MDOF  semidiscrete  structural  models  to  earthquake  excitation.  Along  with 
others,  Chopra  (1995,  page  566)  observes  that  conditionally  stable  procedures  can 
be  used  effectively  for  analysis  of  linear  response  of  large  MDOF  semidiscrete 
structural  models  by  time  history-modal  response  analysis.  This  flexibility  is  due 
to  the  fact  that  only  those  (lower)  modes  that  contribute  to  the  MDOF  model 
response  are  typically  used  in  the  time-history  response  computations  and,  thus, 
control  the  assignment  of  the  largest  time-step  A  t  that  can  be  used  in  the  numerical 
analysis  (Chopra  1995,  pages  574-575,  or  Hughes  1987,  page  493).  Chopra 
(1995,  page  575)  observes  that  when  using  a  direct  solution  of  the  equations  of 
motion  of  a  large,  semidiscrete  MDOF  system  model  or  when  all  modes  are 
included  in  a  time  history-modal  response  analysis,  the  limiting  time-step 
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associated  with  conditionally  stable  algorithms  may  be  prohibitive.  Uncondition¬ 
ally  stable  numerical  methods  are  usually  more  advantageous  for  these  response 
analyses. 


5.4.2  Accuracy  of  response 

All  numerical  step-by-step  procedures  used  in  solving  linear  MDOF  system 
models  must  be  cognizant  of  the  issue  of  accuracy,  whether  the  formulation  used 
is  a  direct  solution  of  the  equation  of  motion  (the  matrix  form  of  Equation  5  in 
Chapter  2)  or  when  using  time  history-modal  response  analysis.  This  numerical 
investigation  of  the  accuracy  of  response  of  SDOF  systems  to  base  excitation  did 
not  include  MDOF  system  models.  Additionally,  it  is  recognized  that  there  are 
many  different  types  of  MDOF  system  response  analysis  formulations  and  eadh  is 
likely  to  have  its  own  unique  characteristics  with  regard  to  the  issue  of  accuracy  of 
computed  response(s).  Some  of  the  factors  contributing  to  these  differences 
include  not  only  the  general  solution  formulation  to  the  equation  of  motion  for  the 
MDOF  system  model,  but  also  the  finite  element  stiffness  and  mass  formulations 
used  in  the  model,  as  well  as  the  method  used  to  incorporate  damping  in  the  analy¬ 
sis.  In  the  interim,  and  until  the  results  from  additional  detailed  numerical  studies 
become  available,  the  following  approximation  is  suggested  to  answer  the  ques¬ 
tion:  Is  the  time-step  in  which  the  ground  acceleration  time-history  is  discretized 
sufficiently  small  to  provide  for  an  adequate  level  of  accuracy  in  the  computed 
dynamic  structural  response(s)  for  the  damped  semidiscrete  MDOF  system  model 
(P  =  0.05)? 

a.  Make  an  initial  selection  of  a  ground  acceleration  time-history  for  the  proj¬ 
ect  and  make  available  information  regarding  the  frequencies  contained 
within  the  acceleration  time-history. 

b.  Develop  a  semidiscrete  or  discrete  (finite  element)  model  of  the  hydraulic 
structure  and  conduct  a  modal  analysis  to  identify  the  first  J  modes  with  a 
significant  contribution  to  system  response  for  the  MDOF  model  compris¬ 
ing  N  modes  (with  J  <  N).  For  many  types  of  hydraulic  structures  the 
dynamic  response  is  often  dominated  by  the  first  few  modes  (i.e., 

J  «<  N).  A  sufficient  number  of  modes  have  been  identified  if  the  sum  of 
the  modal  mass  participation  factors  is  greater  than  or  equal  to  90  percent. 
Alternatively,  simplified  procedures  may  be  used  to  approximately  com¬ 
pute  these  modal  periods  (when  these  procedures  are  available  for  the  type 
of  hydraulic  structure  being  analyzed).  Use  of  simplified  procedures  is 
entirely  appropriate  in  this  exercise  to  compute  the  natural  periods  of  these 
J  modes.  For  example,  Fenves  and  Chopra  (1986)  describe  a  simplified 
procedure  to  compute  the  natural  period  of  a  concrete  gravity  dam  section. 
Experience  has  shown  that  dynamic  response  of  concrete  gravity  dams  is 
dominated  by  first  mode  response  (J  =  1).  A  simplified  procedure  for  com¬ 
puting  the  first  two  modal  periods  (J  =  2)  of  intake  towers  (developed  by 
the  third  and  first  authors  of  this  report)  is  described  in  Appendix  B  of 
Engineer  Circular  1110-2-285  (Headquarters,  U.S.  Army  Corps  of 
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Engineers,  1995).  Experience  has  shown  that  dynamic  response  of  intake 
towers  is  dominated  by  the  first  two  modes. 

c.  Given  the  frequency  content  of  the  ground  motion  fg  of  interest  and  the  At 
to  which  the  acceleration  time-history  is  discretized,  use  the  error  sum¬ 
maries  given  in  Tables  1  through  6  and  in  Figures  13  through  20  to  assess 
the  accuracy  of  computed  response(s)  for  each  of  the  J  modes.  Since  the 
data  given  are  only  for  two  natural  periods  (i.e.,  T0  =  0.25  sec  and  Ta  = 

0.5  sec),  interpolation/extrapolation  may  be  required. 

d.  For  hydraulic  structures  whose  dynamic  response  is  dominated  by  first 
mode  response,  the  error  computed  in  step  3  (c)  is  the  approximation  for 
the  accuracy  of  computed  response  for  the  time-step  At  for  which  the 
ground  acceleration  time-history  is  discretized.  An  assessment  of  the  ade¬ 
quacy  of  the  error  in  computed  response  is  then  made.  If  deemed  accept¬ 
able,  step  3  is  repeated  for  all  other  response  parameters  of  interest.  If  any 
response  parameter  error  is  deemed  unacceptable,  another  ground  accel¬ 
eration  time-history  possessing  a  finer  time-step  should  be  obtained. 

Recall  that  the  time-step  At  for  ground  acceleration  time-histories  is  usu¬ 
ally  equal  to  0.02,  0.01,  or  0.005  sec.  The  process  should  be  repeated 
using  this  new  acceleration  time-history. 

e.  For  hydraulic  structures  whose  response  is  dominated  by  contributions  of 
more  than  one  mode,  the  results  from  the  following  two  simplified  proce¬ 
dures  should  be  considered:  (1)  approximate  the  total  error  in  the  specified 
response  parameter  as  the  largest  of  the  response  errors  of  the  modes;  or 
(2)  approximate  the  total  error  as  the  weighted  sum  of  the  error  of  each 
mode.  The  value  assigned  to  each  weighting  factor  is  intended  to  account 
for  contribution  of  that  mode  to  the  total  response.  Values  of  the  mass 
modal  participation  factors  are  expected  to  be  useful  data  in  this 
evaluation. 


5.4.3  Numerical  damping  of  high-frequency  response 

Numerical  damping  in  a  step-by-step  analysis  of  a  semidiscrete  MDOF  struc¬ 
tural  model  is  advantageous  because  it  filters  out  response  contributions  from 
high-frequency  modes  that  result  from  the  numerical  structural  model  idealization 
and  not  an  actual  property  of  the  structure.  Chopra  (1995,  page  576)  observes 
that  ‘Wilson’s  method  provides  for  numerical  damping  in  modes  with  period  Tn 
such  that  At/Tn  >1.0;  other  methods  are  also  available.”  Refer  to  Chopra  ( 1 995, 
pages  575-576)  or  to  Hughes  (1987,  pages  498-499)  for  further  details  regarding 
this  topic. 


72 


Chapter  5  Results  and  Conclusions 


5.4.4  Nonlinear  analysis 


Classical  modal  analysis  is  not  used  for  time-history  response  analysis  of  non¬ 
linear  structural  systems  because  of  coupling  between  modal  equations  (Chopra 
1995,  page  574).  Unconditionally  stable  procedures  are  generally  necessary  for 
nonlinear  response  analysis  of  such  systems  (Chopra  1995,  page  566). 
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Appendix  A 

Exact  Solution  to  SDOF  System 
Sine  Wave  Base  Excitation 


The  equation  of  motion  for  a  single-degree-of-freedom  system  subjected  to  a 
base  excitation  is 

x(t)  +  2(3cni(0  +  w2x(t)  =  -xp0XSni(t)  (Al) 

where 


Vund(0  ^pgasmat  (A2) 

and  pga  is  the  peak  groimd  acceleration  or  amplitude  of  sinusoidal  base  excitation. 
The  solution  to  this  differential  equation  is 

x(t)  =  e  ^wt(A  cos  (x>Dt  +BsmwDt) 

(1  -/-2)sinG)t  -2$r coscat 
(1  -  r2)2  +  (2pr)2 

where 


0) 


and 
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Al 


(e.g.,  Clough  and  Potmen  1993).  The  values  assigned  to  Tg  and pga.  the  peak 
ground  acceleration,  in  this  study  are  reported  in  Figure  7,  Chapter  4.  Taking  the 
time  derivative  of  Equation  A3  gives 


x(t)  =  e  [(Bu>d  -  AcoP)  cos  u>Dt  -  (5g>P  +  Aa>D)  sincu^t] 


Pga  r 


G) 


(1  -r2)  cos  cot  +  2 Pr  sin  cut 


(r2  -  l)2  +  (2Pr)2 


(A4) 


The  constants  A  and  B  can  be  determined  by  evaluating  Equations  A3  and  A4  for 
the  boundary  conditions  x(t  =  0)  =  0,  and  x(t  =  0)  =  0: 


A  =  Pga _ 2P r _ 

co2  (1  -  r2)1  +  (2Pr)2 


(A5) 


and 


B  =  — — 

(On 


v4coP  + 


pgar  1  -  r2 

(o  (2rP)2  +  (r2-  l)2 


(A6) 


The  relative  acceleration  is  determined  by  taking  the  time  derivative  of  the  relative 
velocity: 


x(t)  =  e- 


A  co2  P2  -  2  Bud  co  P  -  Ao)d  I  cos  0)Dt 


+  B co2p2  +  2^4  (x>D co P  -  Bw>d  sincoD t 


(A7) 


-  pga  rl 


(1  -  r2)  since?  -  2Pr  sin  cot 


(r2-l)2  +  (2Pr)2 


The  total  acceleration,  xtotal  (?),  is  simply  the  sum  of  the  relative  acceleration  plus 
the  ground  acceleration 


WO  =  Xground(0  +  *(0 


=  -[2pcux(r)  +  co2x(r)] 


(AS) 


A2 
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Appendix  B 
Fourier  Series 


Appendix  B 


Fourier  showed  that  any  periodic  function  may  be  expressed  as  the  summation 
of  an  infinite  number  of  sine  and  cosine  terms  (e.g.,  Paz  1985).  Such  a  sum  is 
known  as  a  Fourier  series: 


P(t)  =  aQ  +  ax  cos  wt  +  a2cos2(ot  +  ...  +  ancosmot 

+  b.sincot  +  &,sin2a>t  +  ...  +  b  sinncut 

1  Z  n 


(Bl) 


or 


P(t)  =  aQ  +  ( an  cos  +  bn  sin  nwt)  (B2) 

n  =  1 

where  cb=  2 n/Tg  is  the  circular  frequency,  and  Tg  is  the  period  of  the  base  excita¬ 
tion.  The  evaluation  of  the  coefficients  a0,  a„,  and  bn  is  given  by 


T_ 


f  P(t)dt 
h 


2  _ 

an  =  —  J  P(t)  cos  noitdt 

Ts  h 


2  6  _ 

bn  =  —  J  P(t)  sin  mutdt 

Ts  f, 


Fourier  Series 


(B3) 


(B4) 


(B5) 


r 


B1 


where  tx  in  the  limits  of  the  integrals  may  be  any  value  of  time,  but  is  usually  equal 
-TJ2  or  zero.  The  constant  a0  represents  the  average  of  the  periodic  function  Pit). 

Representing  the  forcing  function  by  a  piecewise  linear  function,  the  Fourier 
coefficients  are  the  summation  of  the  integrals  evaluated  for  each  linear  segment  of 
the  forcing  function: 


(B6) 


N 


an  =  —  f  P(t)  cos  matdt 
T„  ;=i  y 


(B7) 


S  t. 


bn  =  £  f  p(0  sm  nutdt 

T*  -  L 


(B8) 


where  N  =  the  number  of  segments  in  the  piecewise  forcing  function.  The  forcing 
function  in  any  interval  <  t  <  t,  is  expressed  by 


A  t 

where 


(B9) 


AF,  =  Pit)  -  Pit 


(BIO) 


Substituting  Equation  B9  into  Equation  B7  and  performing  the  integration  gives: 


1  a 

=  J-E 


A  t 


Pit)  * Pit ,,!> 


(Bll) 


B2 
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mo 


A P. 


(sinncnr  -  sinrzcor^) 
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+  — — - — [(cos  mot -  cos  mo thl) 
n2(o2At 


(B12) 


+  mo  (/;sin  moti  -  thl  sin went  _, )]  > 
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—  E  \— 

TSU 
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AP 


(cosmott  -  cosmothl) 


+  - - — r(sin«cnr  -  sin«m/.  ,) 

-> — 2  .  ' 
n2(0  At 


(B13) 


+  mo  (tj cos  moti  -  thl  cos  mothl)]  f 


The  response  of  the  SDOF  system  to  a  periodic  force  represented  by  a  Fourier 
series  is  the  superposition  of  the  response  to  each  component  of  the  series.  When 
the  transient  is  omitted,  the  total  steady  state  response  of  a  damped  SDOF  system 
may  be  expressed: 


x.= 


k 


(l-'-„2)2*(2'-„P)2 


sin  motj 
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a  j 

n ' 
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of  y  and  p  that  correspond  to  the  linear  acceleration  method),  the  Wilson  0  Method,  the  Central  Difference  Method,  the 
4th-order  Runge-Kutta  method,  Duhamel's  integral  solved  in  a  piecewise  exact  fashion,  and  the  piecewise  exact  method 
applied  directly. 

Much  of  the  current  guidance  for  selecting  the  time-step  used  in  computing  dynamic  response  to  ground  acceleration  is 
based  on  studies  of  the  accuracy  of  numerical  methods  for  computing  free  vibration  response  of  undamped  SDOF  systems. 
The  information  from  the  free  vibration  studies  is  often  combined  with  useful,  but  qualitative ,  reference  to  the  frequency 
characteristics  of  the  forcing  function.  The  time-step  criteria  for  a  specified  level  of  accuracy  are  often  expressed  as  a  frac¬ 
tion  of  the  natural  (undamped)  period  of  the  SDOF  system.  This  report  quantifies  how  the  accuracy  of  the  six  numerical 
step-by-step  procedures  is  affected  by  the  time-step  and  the  ground  motion  frequency  characteristics. 


